Repository: BQSKit/bqskit (language: python)

Issue:
Substitute



Discussion / hints:

Fix patch:
diff --git a/.github/workflows/pre-commit.yml b/.github/workflows/pre-commit.yml
index b55603d1c..3fa4ed41d 100644
--- a/.github/workflows/pre-commit.yml
+++ b/.github/workflows/pre-commit.yml
@@ -12,7 +12,7 @@ jobs:
     - uses: actions/checkout@v1
     - uses: actions/setup-python@v1
       with:
-        python-version: '3.9'
+        python-version: '3.10'
     - name: set PY
       run: echo "name=PY::$(python -c 'import hashlib, sys;print(hashlib.sha256(sys.version.encode()+sys.executable.encode()).hexdigest())')" >> $GITHUB_ENV
     - uses: actions/cache@v1
diff --git a/.github/workflows/tox.yml b/.github/workflows/tox.yml
index d65d6ae59..509ad9de6 100644
--- a/.github/workflows/tox.yml
+++ b/.github/workflows/tox.yml
@@ -9,7 +9,7 @@ jobs:
   tox:
     strategy:
       matrix:
-        python-version: ['3.7', '3.8', '3.9']
+        python-version: ['3.7', '3.8', '3.9', '3.10']
         os: [macOS-latest, windows-latest, ubuntu-latest]
 
     runs-on: ${{ matrix.os }}
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index bcb86eb3a..e1db1354d 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -1,6 +1,6 @@
 repos:
 -   repo: https://github.com/pre-commit/pre-commit-hooks
-    rev: v4.0.1
+    rev: v4.1.0
     hooks:
     -   id: trailing-whitespace
     -   id: end-of-file-fixer
@@ -30,7 +30,7 @@ repos:
             - --wrap-summaries=80
             - --wrap-descriptions=80
 -   repo: https://github.com/pre-commit/mirrors-autopep8
-    rev: v1.5.7
+    rev: v1.6.0
     hooks:
     -   id: autopep8
         args:
@@ -38,14 +38,14 @@ repos:
             - -a
             - --max-line-length=80
 -   repo: https://github.com/asottile/pyupgrade
-    rev: v2.29.0
+    rev: v2.31.0
     hooks:
     -   id: pyupgrade
         args:
             - --py36-plus
             - --keep-percent-format
 -   repo: https://github.com/asottile/reorder_python_imports
-    rev: v2.6.0
+    rev: v2.7.1
     hooks:
     -   id: reorder-python-imports
         args:
@@ -53,13 +53,13 @@ repos:
             - from __future__ import annotations
             - --py3-plus
 -   repo: https://github.com/asottile/add-trailing-comma
-    rev: v2.1.0
+    rev: v2.2.1
     hooks:
     -   id: add-trailing-comma
         args:
             - --py36-plus
 -   repo: https://github.com/asottile/setup-cfg-fmt
-    rev: v1.18.0
+    rev: v1.20.0
     hooks:
     -   id: setup-cfg-fmt
 - repo: https://github.com/myint/autoflake
@@ -68,11 +68,11 @@ repos:
     - id: autoflake
       args:
         - --in-place
--   repo: https://github.com/pre-commit/mirrors-mypy
-    rev: v0.910-1
-    hooks:
-    - id: mypy
-      additional_dependencies: ["numpy>=1.20,<1.21"]
+# -   repo: https://github.com/pre-commit/mirrors-mypy
+#     rev: v0.931
+#     hooks:
+#     - id: mypy
+#       additional_dependencies: ["numpy>=1.21"]
 -   repo: https://gitlab.com/pycqa/flake8
     rev: 3.9.2
     hooks:
diff --git a/bqskit/compiler/compiler.py b/bqskit/compiler/compiler.py
index 19e968cd0..8dd4d28cb 100644
--- a/bqskit/compiler/compiler.py
+++ b/bqskit/compiler/compiler.py
@@ -38,6 +38,9 @@ def __init__(self, *args: Any, **kwargs: Any) -> None:
             All arguments are passed directly to Dask. You can use
             these to connect to and configure a Dask cluster.
         """
+        if 'silence_logs' not in kwargs:
+            kwargs['silence_logs'] = logging.getLogger('bqskit').level
+
         self.client = Client(*args, **kwargs)
         self.tasks: dict[uuid.UUID, Future] = {}
         _logger.info('Started compiler process.')
@@ -60,12 +63,12 @@ def close(self) -> None:
             self.client.close()
             self.tasks = {}
             _logger.info('Stopped compiler process.')
-        except AttributeError:
+        except (AttributeError, TypeError):
             pass
 
     def submit(self, task: CompilationTask) -> None:
         """Submit a CompilationTask to the Compiler."""
-        executor = Executor(task)
+        executor = self.client.scatter(Executor(task))
         future = self.client.submit(Executor.run, executor, pure=False)
         self.tasks[task.task_id] = future
         _logger.info('Submitted task: %s' % task.task_id)
@@ -76,7 +79,7 @@ def status(self, task: CompilationTask) -> str:
 
     def result(self, task: CompilationTask) -> Circuit:
         """Block until the CompilationTask is finished, return its result."""
-        circ = self.tasks[task.task_id].result()
+        circ = self.tasks[task.task_id].result()[0]
         return circ
 
     def cancel(self, task: CompilationTask) -> None:
@@ -90,3 +93,9 @@ def compile(self, task: CompilationTask) -> Circuit:
         self.submit(task)
         result = self.result(task)
         return result
+
+    def analyze(self, task: CompilationTask, key: str) -> Any:
+        """Gather the value associated with `key` in the task's data."""
+        if task.task_id not in self.tasks:
+            self.submit(task)
+        return self.tasks[task.task_id].result()[1][key]
diff --git a/bqskit/compiler/executor.py b/bqskit/compiler/executor.py
index bd86e67af..4a1517f7c 100644
--- a/bqskit/compiler/executor.py
+++ b/bqskit/compiler/executor.py
@@ -31,10 +31,10 @@ def __init__(self, task: CompilationTask) -> None:
         self.data: dict[str, Any] = {'executor': self}
         self.done = False
 
-    def run(self) -> Circuit:
+    def run(self) -> tuple[Circuit, dict[str, Any]]:
         """Execute the task."""
         with threadpool_limits(limits=1):
             for pass_obj in self.passes:
                 pass_obj.run(self.circuit, self.data)
         self.done = True
-        return self.circuit
+        return self.circuit, self.data
diff --git a/bqskit/compiler/machine.py b/bqskit/compiler/machine.py
index 3c6e9bf4d..a01685863 100644
--- a/bqskit/compiler/machine.py
+++ b/bqskit/compiler/machine.py
@@ -94,7 +94,7 @@ def get_locations(self, block_size: int) -> list[CircuitLocation]:
             )
 
         if block_size <= 0:
-            raise ValueError('Nonpositive block_size; got {block_size}.')
+            raise ValueError(f'Nonpositive block_size; got {block_size}.')
 
         locations: set[CircuitLocation] = set()
 
diff --git a/bqskit/compiler/task.py b/bqskit/compiler/task.py
index abaaee2b9..93600eecc 100644
--- a/bqskit/compiler/task.py
+++ b/bqskit/compiler/task.py
@@ -7,16 +7,17 @@
 
 from bqskit.compiler.basepass import BasePass
 from bqskit.ir.circuit import Circuit
-from bqskit.ir.gates.constant.cx import CNOTGate
+from bqskit.ir.gates import CircuitGate
+from bqskit.ir.gates import CNOTGate
+from bqskit.ir.operation import Operation
 from bqskit.passes.control import ForEachBlockPass
 from bqskit.passes.control.predicates.count import GateCountPredicate
 from bqskit.passes.control.whileloop import WhileLoopPass
 from bqskit.passes.partitioning.cluster import ClusteringPartitioner
 from bqskit.passes.partitioning.quick import QuickPartitioner
 from bqskit.passes.processing import ScanningGateRemovalPass
-from bqskit.passes.processing import WindowOptimizationPass
 from bqskit.passes.synthesis import QFASTDecompositionPass
-from bqskit.passes.synthesis.leap import OptimizedLEAPPass
+from bqskit.passes.synthesis.leap import LEAPSynthesisPass
 from bqskit.passes.util import UnfoldPass
 from bqskit.qis.unitary.unitarymatrix import UnitaryLike
 
@@ -57,15 +58,19 @@ def synthesize(utry: UnitaryLike) -> CompilationTask:
             _logger.warning('Synthesis input size is very large.')
 
         inner_seq = [
-            OptimizedLEAPPass(),
-            WindowOptimizationPass(),
+            LEAPSynthesisPass(),
             ScanningGateRemovalPass(),
         ]
 
         passes: list[BasePass] = []
         if num_qudits >= 5:
             passes.append(QFASTDecompositionPass())
-            passes.append(ForEachBlockPass(inner_seq))
+            passes.append(
+                ForEachBlockPass(
+                    inner_seq,
+                    replace_filter=less_2q_gates,
+                ),
+            )
             passes.append(UnfoldPass())
         else:
             passes.extend(inner_seq)
@@ -81,23 +86,32 @@ def optimize(circuit: Circuit) -> CompilationTask:
             return CompilationTask.synthesize(circuit.get_unitary())
 
         inner_seq = [
-            OptimizedLEAPPass(),
+            LEAPSynthesisPass(),
             ScanningGateRemovalPass(),
         ]
 
         passes: list[BasePass] = []
         passes.append(QuickPartitioner(3))
-        passes.append(ForEachBlockPass(inner_seq))
+        passes.append(ForEachBlockPass(inner_seq, replace_filter=less_2q_gates))
         passes.append(UnfoldPass())
 
         iterative_reopt = WhileLoopPass(
             GateCountPredicate(CNOTGate()),
             [
                 ClusteringPartitioner(3, 4),
-                ForEachBlockPass(inner_seq),
+                ForEachBlockPass(inner_seq, replace_filter=less_2q_gates),
                 UnfoldPass(),
             ],
         )
 
         passes.append(iterative_reopt)
         return CompilationTask(circuit, passes)
+
+
+def less_2q_gates(circuit: Circuit, op: Operation) -> bool:
+    """Replace `circuit' with `op` if has less 2 qubit gates."""
+    if not isinstance(op, CircuitGate):
+        return True
+    og_num_2q_gate = len([op for op in op._circuit if op.num_qudits >= 2])
+    new_num_2q_gate = len([op for op in circuit if op.num_qudits >= 2])
+    return new_num_2q_gate > og_num_2q_gate
diff --git a/bqskit/exec/__init__.py b/bqskit/exec/__init__.py
new file mode 100644
index 000000000..4f7c0497e
--- /dev/null
+++ b/bqskit/exec/__init__.py
@@ -0,0 +1,47 @@
+"""
+=============================================================
+Circuit Execution Module (:mod:`bqskit.exec`)
+=============================================================
+
+.. currentmodule:: bqskit.exec
+
+BQSKit supports a variety of ways to execute circuits through
+CircuitRunner.
+
+.. rubric:: Core Classes
+
+.. autosummary::
+    :toctree: autogen
+    :recursive:
+
+    CircuitRunner
+    RunnerResults
+
+.. rubric:: CircuitRunners
+
+.. autosummary::
+    :toctree: autogen
+    :recursive:
+
+    QuestRunner
+    IBMQRunner
+    SimulationRunner
+"""
+from __future__ import annotations
+
+__all__ = [
+    'CircuitRunner',
+    'RunnerResults',
+    'QuestRunner',
+    'SimulationRunner',
+]
+
+from bqskit.exec.results import RunnerResults
+from bqskit.exec.runner import CircuitRunner
+try:
+    from bqskit.exec.runners.ibmq import IBMQRunner
+    __all__ += ['IBMQRunner']
+except ImportError:
+    pass
+from bqskit.exec.runners.quest import QuestRunner
+from bqskit.exec.runners.sim import SimulationRunner
diff --git a/bqskit/exec/results.py b/bqskit/exec/results.py
new file mode 100644
index 000000000..e64a33384
--- /dev/null
+++ b/bqskit/exec/results.py
@@ -0,0 +1,27 @@
+from __future__ import annotations
+
+from typing import Sequence
+
+import numpy as np
+import numpy.typing as npt
+
+
+class RunnerResults:
+    """The result of running a Circuit through a CircuitRunner."""
+
+    def __init__(
+        self,
+        num_qudits: int,
+        radixes: Sequence[int],
+        probs: Sequence[float],
+    ) -> None:
+        """Construct a RunnerResults object."""
+        self.num_qudits = num_qudits
+        self.radixes = radixes
+        self.probs = np.array(probs)
+
+    def get_counts(self, shots: int) -> npt.NDArray[np.int64]:
+        return np.asarray(np.multiply(shots, self.probs), np.int64)
+
+    def __str__(self) -> str:
+        pass
diff --git a/bqskit/exec/runner.py b/bqskit/exec/runner.py
new file mode 100644
index 000000000..2d811a9e4
--- /dev/null
+++ b/bqskit/exec/runner.py
@@ -0,0 +1,23 @@
+from __future__ import annotations
+
+from abc import ABC
+from abc import abstractmethod
+
+from bqskit.exec.results import RunnerResults
+from bqskit.ir.circuit import Circuit
+
+
+class CircuitRunner(ABC):
+    """A CircuitRunner is responsible for executing a quantum circuit."""
+
+    @abstractmethod
+    def run(self, circuit: Circuit) -> RunnerResults:
+        """
+        Execute the circuit and return results.
+
+        Args:
+            circuit (Circuit): The circuit to run.
+
+        Returns:
+            (RunnerResults): The results from executing this circuit.
+        """
diff --git a/bqskit/exec/runners/__init__.py b/bqskit/exec/runners/__init__.py
new file mode 100644
index 000000000..1377ef5d4
--- /dev/null
+++ b/bqskit/exec/runners/__init__.py
@@ -0,0 +1,13 @@
+"""This package contains several CircuitRunner implementations."""
+from __future__ import annotations
+
+__all__ = ['QuestRunner', 'SimulationRunner']
+
+try:
+    from bqskit.exec.runners.ibmq import IBMQRunner
+    __all__ += ['IBMQRunner']
+except ImportError:
+    pass
+
+from bqskit.exec.runners.quest import QuestRunner
+from bqskit.exec.runners.sim import SimulationRunner
diff --git a/bqskit/exec/runners/ibmq.py b/bqskit/exec/runners/ibmq.py
new file mode 100644
index 000000000..c9735da60
--- /dev/null
+++ b/bqskit/exec/runners/ibmq.py
@@ -0,0 +1,34 @@
+from __future__ import annotations
+
+from qiskit import QuantumCircuit
+from qiskit.providers.ibmq import IBMQBackend
+
+from bqskit.exec.results import RunnerResults
+from bqskit.exec.runner import CircuitRunner
+from bqskit.ir.circuit import Circuit
+
+
+class IBMQRunner(CircuitRunner):
+    """Simulate a circuit."""
+
+    def __init__(self, backend: IBMQBackend) -> None:
+        """Setup an IBMQRunner to execute circuits on `backend`."""
+        self.backend = backend
+
+    def run(self, circuit: Circuit) -> RunnerResults:
+        """Execute the circuit, see CircuitRunner.run for more info."""
+        # 1. Check circuit and self.backend are compatible
+        # TODO
+
+        # 2. Convert to Qiskit IR
+        qiskit_circ = QuantumCircuit.from_qasm_str(circuit.to('qasm'))
+        qiskit_circ.measure_all()
+
+        # 3. Run circuit
+        result = self.backend.run(qiskit_circ).result()
+        shots = result.results[0].shots
+        probs = [0.0 for i in range(2 ** circuit.num_qudits)]
+        for bit_str, count in result.get_counts().items():
+            probs[int(bit_str, 2)] = count / shots
+
+        return RunnerResults(circuit.num_qudits, circuit.radixes, probs)
diff --git a/bqskit/exec/runners/quest.py b/bqskit/exec/runners/quest.py
new file mode 100644
index 000000000..8a34ba32b
--- /dev/null
+++ b/bqskit/exec/runners/quest.py
@@ -0,0 +1,280 @@
+"""This module implements the QuestRunner CircuitRunner."""
+from __future__ import annotations
+
+import logging
+import time
+from typing import Any
+
+import numpy as np
+import numpy.typing as npt
+from scipy.optimize import dual_annealing
+
+from bqskit.compiler import CompilationTask
+from bqskit.compiler import Compiler
+from bqskit.exec.results import RunnerResults
+from bqskit.exec.runner import CircuitRunner
+from bqskit.ir.circuit import Circuit
+from bqskit.ir.gates import CircuitGate
+from bqskit.ir.gates import CNOTGate
+from bqskit.ir.operation import Operation
+from bqskit.ir.point import CircuitPoint
+from bqskit.passes.control import ForEachBlockPass
+from bqskit.passes.partitioning import QuickPartitioner
+from bqskit.passes.synthesis import LEAPSynthesisPass
+
+_logger = logging.getLogger(__name__)
+
+
+class QuestRunner(CircuitRunner):
+    """
+    Simulate a circuit using the Quest algorithm.
+
+    References:
+        Patel, Tirthak, et al. "Robust and Resource-Efficient Quantum Circuit
+        Approximation." arXiv preprint arXiv:2108.12714 (2021).
+    """
+
+    def __init__(
+        self,
+        sub_runner: CircuitRunner,
+        block_size: int = 3,
+        compiler: Compiler | None = None,
+        weit: float = 0.5,
+        pert: int = 100,
+        approx_threshold: float | None = None,
+        sample_size: int = 16,
+    ) -> None:
+        """
+        Run the QUEST algorithm using `sub_runner` to execute circuits.
+
+        Args:
+            sub_runner (CircuitRunner): The CircuitRunner to execute the
+                approximated circuits.
+
+            block_size (int): The maximum number of qudits in each
+                block. Defaults to 3.
+
+            compiler (Compiler | None): The compiler object used
+                to partitioning and synthesize the circuit. If left
+                as None, then a standard compiler will be created.
+
+            weit (float): The weight set on number of CNOTs versus
+                approximation dissimilarity. Default 0.5.
+
+            pert (int): The percentile of sample distances used to judge
+                approximation quality.
+
+            approx_threshold (float | None): The maximum allowed error
+                in the approximation. If left as None, then it will be
+                computed based on the number of blocks.
+
+            sample_size (int): The number of approximations to generate.
+        """
+        self.sub_runner = sub_runner
+        self.block_size = block_size
+        self.compiler = compiler if compiler is not None else Compiler()
+        self.weit = weit
+        self.pert = pert
+        self.approx_threshold = approx_threshold
+        self.sample_size = sample_size
+
+    def run(self, circuit: Circuit) -> RunnerResults:
+        """Execute the circuit, see CircuitRunner.run for more info."""
+
+        # 1. Compile the circuit
+        synthesis_pass = LEAPSynthesisPass(store_partial_solutions=True)
+        task = CompilationTask(
+            circuit.copy(), [
+                QuickPartitioner(self.block_size),
+                ForEachBlockPass(synthesis_pass),
+            ],
+        )
+        blocked_circuit = self.compiler.compile(task)
+
+        # 2. Gather partial solutions
+        data = self.compiler.analyze(task, ForEachBlockPass.key)
+        psols, pts = self.parse_data(blocked_circuit, data)
+        # psols: psols[i] = list[[circuit, dist]] -> block i's partial solutions
+        # pts: pts[i] = CircuitPoint -> block i's locations
+
+        # 3. Approximate circuit
+        approx_circuits = self.approximate_circuit(blocked_circuit, psols, pts)
+
+        # 4. Run circuits
+        results = [self.sub_runner.run(c) for c in approx_circuits]
+
+        # 5. Average and return results
+        probs = np.sum(np.array([result.probs for result in results]), axis=0)
+        probs /= self.sample_size
+        return RunnerResults(circuit.num_qudits, circuit.radixes, probs)
+
+    def parse_data(
+        self,
+        blocked_circuit: Circuit,
+        data: dict[Any, Any],
+    ) -> tuple[list[list[tuple[Circuit, float]]], list[CircuitPoint]]:
+        """Parse the data outputed from synthesis."""
+        block_data = data[0]
+
+        psols: list[list[tuple[Circuit, float]]] = []
+        pts: list[CircuitPoint] = []
+
+        for block in block_data:
+            pts.append(block['point'])
+            exact_block = blocked_circuit[pts[-1]].gate._circuit.copy()  # type: ignore  # noqa
+            exact_block.set_params(blocked_circuit[pts[-1]].params)
+            exact_utry = exact_block.get_unitary()
+            psols.append([(exact_block, 0.0)])
+
+            if 'psols' not in block:
+                continue
+
+            for depth, psol_list in block['psols'].items():
+                for psol in psol_list:
+                    dist = psol[0].get_unitary().get_distance_from(exact_utry)
+                    psols[-1].append((psol[0], dist))
+
+        return psols, pts
+
+    def approximate_circuit(
+        self,
+        circuit: Circuit,
+        psols: list[list[tuple[Circuit, float]]],
+        pts: list[CircuitPoint],
+    ) -> list[Circuit]:
+        """Use partial block solutions to compute circuit approximations."""
+        num_blocks = len(pts)
+        approx_threshold = self.approx_threshold or num_blocks / 10
+        bounds = [[0, len(psol_list)] for psol_list in psols]
+
+        psol_configs: list[list[int]] = []
+        approx_circuits: list[Circuit] = []
+        dists = [[psol[1] for psol in psol_list] for psol_list in psols]
+        utrys = [
+            [psol[0].get_unitary() for psol in psol_list]
+            for psol_list in psols
+        ]
+        n_cxs = [
+            [psol[0].count(CNOTGate()) for psol in psol_list]
+            for psol_list in psols
+        ]
+        n_cx = sum(op.gate._circuit.count(CNOTGate()) for op in circuit)  # type: ignore  # noqa
+
+        for _ in range(self.sample_size):
+            x = []
+            f = float('inf')
+            for _ in range(3):
+                result = dual_annealing(
+                    annealing_objective,
+                    args=(
+                        approx_threshold,
+                        n_cx,
+                        psol_configs,
+                        dists,
+                        utrys,
+                        n_cxs,
+                        self.pert,
+                        self.weit,
+                    ),
+                    x0=[0] * num_blocks,
+                    bounds=bounds,
+                    seed=int(time.time()),
+                    initial_temp=5e4,
+                    maxiter=10000,
+                )
+
+                if result.fun < f:
+                    x = result.x
+                    f = result.fun
+
+            blocks = [int(x_e) for x_e in x]
+
+            if blocks in psol_configs:
+                _logger.debug('Generated a repeat approximate circuit.')
+                break
+
+            approx_dist = sum(
+                psols[b][ind][1]
+                for b, ind in enumerate(blocks)
+            )
+
+            if approx_dist > approx_threshold:
+                _logger.debug('Approximate circuit has high distance.')
+                break
+
+            psol_configs.append(blocks)
+            approx_circuit = self.assemble_circuit(circuit, blocks, psols, pts)
+            approx_circuits.append(approx_circuit)
+            _logger.info(
+                'Generated approximate circuit with approximate distance:'
+                f' {approx_dist}.',
+            )
+            _logger.info(
+                f'Approximate circuit has {approx_circuit.count(CNOTGate())}'
+                ' cnots.',
+            )
+
+        return approx_circuits
+
+    def assemble_circuit(
+        self,
+        circuit: Circuit,
+        blocks: list[int],
+        psols: list[list[tuple[Circuit, float]]],
+        pts: list[CircuitPoint],
+    ) -> Circuit:
+        """Assemble a circuit from a list of block indices."""
+        circuit_gates = [
+            CircuitGate(psol_list[b][0])
+            for b, psol_list
+            in zip(blocks, psols)
+        ]
+        locations = [circuit[pt].location for pt in pts]
+        operations = [
+            Operation(cg, loc, cg._circuit.params)
+            for cg, loc
+            in zip(circuit_gates, locations)
+        ]
+        copied_circuit = circuit.copy()
+        copied_circuit.batch_replace(pts, operations)
+        return copied_circuit
+
+
+def annealing_objective(x: npt.NDArray[np.float64], *args: Any) -> float:
+    """Compute the objective function used during approximation generation."""
+    blocks = [int(ind) for ind in x]
+
+    approx_threshold = args[0]
+    upper_limit = args[1]
+    psol_configs = args[2]
+    dists = args[3]
+    utrys = args[4]
+    n_cxs = args[5]
+    p = args[6]
+    w = args[7]
+
+    if blocks in psol_configs:
+        return 1.2
+
+    approx_dist = sum(dists[i][b] for i, b in enumerate(blocks))
+
+    if approx_dist > approx_threshold:
+        return 1.1
+
+    n_cx = w * (sum(n_cxs[i][b] for i, b in enumerate(blocks)) / upper_limit)
+
+    if len(psol_configs) == 0:
+        return n_cx
+
+    distances = [0] * len(psol_configs)
+
+    for index, config in enumerate(psol_configs):
+        distances[index] = sum(  # type: ignore
+            utrys[i][config[i]].get_distance_from(utrys[i][b])
+            < max(dists[i][config[i]], dists[i][b])
+            for i, b in enumerate(blocks)
+        ) / len(blocks)
+
+    b_dv = (1 - w) * np.percentile(distances, p)
+
+    return n_cx + b_dv
diff --git a/bqskit/exec/runners/sim.py b/bqskit/exec/runners/sim.py
new file mode 100644
index 000000000..14c589cec
--- /dev/null
+++ b/bqskit/exec/runners/sim.py
@@ -0,0 +1,17 @@
+from __future__ import annotations
+
+import numpy as np
+
+from bqskit.exec.results import RunnerResults
+from bqskit.exec.runner import CircuitRunner
+from bqskit.ir.circuit import Circuit
+
+
+class SimulationRunner(CircuitRunner):
+    """Simulate a circuit."""
+
+    def run(self, circuit: Circuit) -> RunnerResults:
+        """Execute the circuit, see CircuitRunner.run for more info."""
+        state = np.reshape(circuit.get_unitary()[:, 0], (-1))
+        probs = np.square(np.abs(state), dtype=np.float64)
+        return RunnerResults(circuit.num_qudits, circuit.radixes, probs)  # type: ignore # noqa
diff --git a/bqskit/ir/circuit.py b/bqskit/ir/circuit.py
index f8b4d382f..d16d307b1 100644
--- a/bqskit/ir/circuit.py
+++ b/bqskit/ir/circuit.py
@@ -15,7 +15,7 @@
 from typing import TYPE_CHECKING
 
 import numpy as np
-from distributed import worker_client
+import numpy.typing as npt
 
 from bqskit.ir.gate import Gate
 from bqskit.ir.gates.circuitgate import CircuitGate
@@ -30,8 +30,6 @@
 from bqskit.ir.opt.cost.functions import HilbertSchmidtCost
 from bqskit.ir.opt.instantiater import Instantiater
 from bqskit.ir.opt.instantiaters import instantiater_order
-from bqskit.ir.opt.instantiaters.minimization import Minimization
-from bqskit.ir.opt.instantiaters.qfactor import QFactor
 from bqskit.ir.opt.minimizers.ceres import CeresMinimizer
 from bqskit.ir.point import CircuitPoint
 from bqskit.ir.point import CircuitPointLike
@@ -46,10 +44,13 @@
 from bqskit.qis.unitary.unitarybuilder import UnitaryBuilder
 from bqskit.qis.unitary.unitarymatrix import UnitaryLike
 from bqskit.qis.unitary.unitarymatrix import UnitaryMatrix
+from bqskit.utils.random import seed_random_sources
 from bqskit.utils.typing import is_integer
 from bqskit.utils.typing import is_iterable
 from bqskit.utils.typing import is_sequence_of_int
+from bqskit.utils.typing import is_square_matrix
 from bqskit.utils.typing import is_valid_radixes
+from bqskit.utils.typing import is_vector
 
 if TYPE_CHECKING:
     from bqskit.ir.opt.cost.function import CostFunction
@@ -179,7 +180,7 @@ def num_cycles(self) -> int:
         return len(self._circuit)
 
     @property
-    def params(self) -> np.ndarray:
+    def params(self) -> npt.NDArray[np.float64]:
         """The stored parameters for the circuit."""
         return np.array(sum((list(op.params) for op in self), []))
 
@@ -1227,7 +1228,7 @@ def replace(self, point: CircuitPointLike, op: Operation) -> None:
 
             ValueError: If `point.qudit` is not in `op.location`
         """
-        if point[1] not in op.location:
+        if len(self[point].location.intersection(op.location)) == 0:
             raise ValueError("Point's qudit is not in operation's location.")
 
         self.pop(point)
@@ -2140,13 +2141,13 @@ def get_statevector(self, in_state: StateLike) -> StateVector:
         # TODO: Can be made a lot more efficient.
         return self.get_unitary().get_statevector(in_state)
 
-    def get_grad(self, params: RealVector = []) -> np.ndarray:
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
         """Return the gradient of the circuit."""
         return self.get_unitary_and_grad(params)[1]
 
     def get_unitary_and_grad(
         self, params: RealVector = [],
-    ) -> tuple[UnitaryMatrix, np.ndarray]:
+    ) -> tuple[UnitaryMatrix, npt.NDArray[np.complex128]]:
         """Return the unitary and gradient of the circuit."""
         if len(params) != 0:
             self.check_parameters(params)
@@ -2203,6 +2204,7 @@ def instantiate(
         target: StateLike | UnitaryLike,
         method: str | None = None,
         multistarts: int = 1,
+        seed: int | None = None,
         **kwargs: Any,
     ) -> Circuit:
         """
@@ -2227,6 +2229,10 @@ def instantiate(
             multistarts (int): The number of instantiation jobs to spawn
                 and manage. (Default: 1)
 
+            seed (int | None): The seed for any pseudo-random number generators
+                to use. Note that this is not guaranteed to make this method
+                reproducible.
+
             kwargs (dict[str, Any]): Method specific options, passed
                 directly to method constructor. For more info, see
                 `bqskit.ir.opt.instantiaters`.
@@ -2242,96 +2248,74 @@ def instantiate(
             ValueError: If `target` dimension doesn't match with circuit.
 
             ValueError: If `multistarts` is not a positive integer.
-        """
-        # Check if should run in parallel
-        parallel = 'parallel' in kwargs and kwargs['parallel']
-        if 'parallel' in kwargs:
-            kwargs.pop('parallel')
-
-        # Assign method if unspecified
-        if method is None:
-            error_msg = ''
-            for inst in instantiater_order:
-                if inst.is_capable(self):  # type: ignore
-                    method = inst.get_method_name()  # type: ignore
-                    break
-                else:
-                    report = inst.get_violation_report(self)  # type: ignore
-                    error_msg += report + '\n'
 
-            if method is None:
+            ValueError: If `seed` is not an integer or `None`
+        """
+        # Set seed if specified
+        if seed is not None:
+            if not isinstance(seed, int):
                 raise ValueError(
-                    'No instantiation method works for this circuit.\n%s'
-                    % error_msg,
+                    f'Expected seed to be an integer got {type(seed)}.',
                 )
+            seed_random_sources(seed)
+
+        # Construct instantiater from method
+        instantiater: Instantiater | None = None
+
+        err = ''
+        for inst in instantiater_order:
+            # If method is specified; match it
+            if inst.get_method_name().lower() == method:  # type: ignore
+                if not inst.is_capable(self):  # type: ignore
+                    raise ValueError(
+                        'Circuit cannot be instantiated using the '
+                        f'{method} method.'
+                        f'\n{inst.get_violation_report(self)}',  # type: ignore
+                    )
+                instantiater = inst(**kwargs)
+                break
 
-        # Create instantiater
-        instantiater: Instantiater
-
-        if method.lower() == 'qfactor':
-            instantiater = QFactor(**kwargs)
+            # If method is not specified; find first capable one
+            if method is None:
+                if inst.is_capable(self):  # type: ignore
+                    instantiater = inst(**kwargs)
+                    break
+                err += inst.get_violation_report(self) + '\n'  # type: ignore
 
-        elif method.lower() == 'minimization':
-            instantiater = Minimization(**kwargs)
+        if instantiater is None:
+            if err != '':
+                raise ValueError(f'No capable instantiater.\n{err}')
+            raise ValueError(f'No such instantiatation method {method}.')
 
+        # Check Target
+        if is_square_matrix(target):
+            target = UnitaryMatrix(target)  # type: ignore
+        elif is_vector(target):
+            target = StateVector(target)  # type: ignore
         else:
-            raise ValueError(
-                'No such method %s; expected "qfactor" or "minimization".',
-            )
-
-        if not instantiater.is_capable(self):
-            raise ValueError(
-                'Circuit cannot be instantiated using the %s method.\n %s'
-                % (method, instantiater.get_violation_report(self)),
+            raise TypeError(
+                'Expected either StateVector or UnitaryMatrix'
+                ' for target, got %s.' % type(target),
             )
 
-        # Check Target
-        try:
-            typed_target = UnitaryMatrix(target)  # type: ignore
-        except (ValueError, TypeError):
-            try:
-                typed_target = StateVector(target)  # type: ignore
-            except (ValueError, TypeError) as ex:
-                raise TypeError(
-                    'Expected either StateVector or UnitaryMatrix'
-                    ' for target, got %s.' % type(target),
-                ) from ex
-
-        if typed_target.dim != self.dim:
+        if target.dim != self.dim:
             raise ValueError('Target dimension mismatch with circuit.')
 
         # Generate starting points
-        starts = instantiater.gen_starting_points(
-            multistarts,
-            self,
-            typed_target,
-        )
+        starts = instantiater.gen_starting_points(multistarts, self, target)
 
         # Instantiate the circuit
-        params = []
-        if multistarts > 1 and parallel:
-            with worker_client() as client:
-                futures = [
-                    client.submit(
-                        instantiater.instantiate,
-                        self,
-                        typed_target,
-                        start,
-                    )
-                    for start in starts
-                ]
-                params = [future.result() for future in futures]
+        params_list = []
+        for start in starts:
+            params_list.append(instantiater.instantiate(self, target, start))
 
-        else:
-            for start in starts:
-                params.append(
-                    instantiater.instantiate(
-                        self, typed_target, start,
-                    ),
-                )
+        # Return best result
+        if multistarts == 1:
+            self.set_params(params_list[0])
+            return self
 
-        cost_fn = HilbertSchmidtCost(self, typed_target)
-        self.set_params(sorted(params, key=lambda x: cost_fn(x))[0])
+        cost_fn = HilbertSchmidtCost(self, target)
+        self.set_params(sorted(params_list, key=lambda x: cost_fn(x))[0])
         return self
 
     def minimize(self, cost: CostFunction, **kwargs: Any) -> None:
@@ -2666,6 +2650,10 @@ def save(self, filename: str) -> None:
         with open(filename, 'w') as f:
             f.write(language.encode(self))
 
+    def to(self, type: str) -> str:
+        """Convert circuit to language."""
+        return get_language(type).encode(self)
+
     @staticmethod
     def from_file(filename: str) -> Circuit:
         """Restore a circuit from a file."""
diff --git a/bqskit/ir/gates/circuitgate.py b/bqskit/ir/gates/circuitgate.py
index 952775efb..5bbc747b3 100644
--- a/bqskit/ir/gates/circuitgate.py
+++ b/bqskit/ir/gates/circuitgate.py
@@ -4,6 +4,7 @@
 from typing import TYPE_CHECKING
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.qis.unitary.unitary import RealVector
 from bqskit.qis.unitary.unitarymatrix import UnitaryMatrix
@@ -44,7 +45,7 @@ def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
         """Return the unitary for this gate, see :class:`Unitary` for more."""
         return self._circuit.get_unitary(params)
 
-    def get_grad(self, params: RealVector = []) -> np.ndarray:
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
         """
         Return the gradient for this gate.
 
@@ -55,7 +56,7 @@ def get_grad(self, params: RealVector = []) -> np.ndarray:
     def get_unitary_and_grad(
         self,
         params: RealVector = [],
-    ) -> tuple[UnitaryMatrix, np.ndarray]:
+    ) -> tuple[UnitaryMatrix, npt.NDArray[np.complex128]]:
         """
         Return the unitary and gradient for this gate.
 
diff --git a/bqskit/ir/gates/composed/controlled.py b/bqskit/ir/gates/composed/controlled.py
index 9758fec15..0d023635a 100644
--- a/bqskit/ir/gates/composed/controlled.py
+++ b/bqskit/ir/gates/composed/controlled.py
@@ -2,6 +2,7 @@
 from __future__ import annotations
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.ir.gate import Gate
 from bqskit.ir.gates.composedgate import ComposedGate
@@ -61,9 +62,9 @@ def __init__(
         self._name = '%d-Controlled(%s)' % (num_controls, gate.name)
         self._num_params = gate.num_params
 
-        self.Ic = np.identity(2 ** num_controls)
-        self.It = np.identity(gate.dim)
-        self.OneProj = np.zeros(self.Ic.shape)
+        self.Ic = np.identity(2 ** num_controls, dtype=np.complex128)
+        self.It = np.identity(gate.dim, dtype=np.complex128)
+        self.OneProj = np.zeros(self.Ic.shape, dtype=np.complex128)
         self.OneProj[-1, -1] = 1
         self.left = np.kron((self.Ic - self.OneProj), self.It)
 
@@ -82,7 +83,7 @@ def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
         right = np.kron(self.OneProj, U)
         return UnitaryMatrix(self.left + right, self.radixes)
 
-    def get_grad(self, params: RealVector = []) -> np.ndarray:
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
         """
         Return the gradient for this gate.
 
@@ -98,7 +99,7 @@ def get_grad(self, params: RealVector = []) -> np.ndarray:
     def get_unitary_and_grad(
         self,
         params: RealVector = [],
-    ) -> tuple[UnitaryMatrix, np.ndarray]:
+    ) -> tuple[UnitaryMatrix, npt.NDArray[np.complex128]]:
         """
         Return the unitary and gradient for this gate.
 
@@ -113,3 +114,13 @@ def get_unitary_and_grad(
         utry = UnitaryMatrix(self.left + right, self.radixes)
         grads = np.kron(self.OneProj, G)
         return utry, grads
+
+    def __eq__(self, other: object) -> bool:
+        return (
+            isinstance(other, ControlledGate)
+            and self.gate == other.gate
+            and self.num_controls == other.num_controls
+        )
+
+    def __hash__(self) -> int:
+        return hash((self.gate, self.num_controls))
diff --git a/bqskit/ir/gates/composed/daggergate.py b/bqskit/ir/gates/composed/daggergate.py
index 7c747cfa1..c38ddcc90 100644
--- a/bqskit/ir/gates/composed/daggergate.py
+++ b/bqskit/ir/gates/composed/daggergate.py
@@ -2,6 +2,7 @@
 from __future__ import annotations
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.ir.gate import Gate
 from bqskit.ir.gates.composedgate import ComposedGate
@@ -54,7 +55,7 @@ def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
 
         return self.gate.get_unitary(params).dagger
 
-    def get_grad(self, params: RealVector = []) -> np.ndarray:
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
         """
         Return the gradient for this gate.
 
@@ -73,7 +74,7 @@ def get_grad(self, params: RealVector = []) -> np.ndarray:
     def get_unitary_and_grad(
         self,
         params: RealVector = [],
-    ) -> tuple[UnitaryMatrix, np.ndarray]:
+    ) -> tuple[UnitaryMatrix, npt.NDArray[np.complex128]]:
         """
         Return the unitary and gradient for this gate.
 
@@ -85,7 +86,7 @@ def get_unitary_and_grad(
         utry, grads = self.gate.get_unitary_and_grad(params)  # type: ignore
         return utry.dagger, np.transpose(grads.conj(), (0, 2, 1))
 
-    def optimize(self, env_matrix: np.ndarray) -> list[float]:
+    def optimize(self, env_matrix: npt.NDArray[np.complex128]) -> list[float]:
         """
         Return the optimal parameters with respect to an environment matrix.
 
@@ -95,3 +96,12 @@ def optimize(self, env_matrix: np.ndarray) -> list[float]:
             return []
         self.check_env_matrix(env_matrix)
         return self.gate.optimize(env_matrix.conj().T)  # type: ignore
+
+    def __eq__(self, other: object) -> bool:
+        return (
+            isinstance(other, DaggerGate)
+            and self.gate == other.gate
+        )
+
+    def __hash__(self) -> int:
+        return hash(self.gate)
diff --git a/bqskit/ir/gates/composed/frozenparam.py b/bqskit/ir/gates/composed/frozenparam.py
index 4488ccf57..53dd7a8a6 100644
--- a/bqskit/ir/gates/composed/frozenparam.py
+++ b/bqskit/ir/gates/composed/frozenparam.py
@@ -2,6 +2,7 @@
 from __future__ import annotations
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.ir.gate import Gate
 from bqskit.ir.gates.composedgate import ComposedGate
@@ -92,7 +93,7 @@ def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
         """Return the unitary for this gate, see :class:`Unitary` for more."""
         return self.gate.get_unitary(self.get_full_params(params))
 
-    def get_grad(self, params: RealVector = []) -> np.ndarray:
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
         """
         Return the gradient for this gate.
 
@@ -104,7 +105,7 @@ def get_grad(self, params: RealVector = []) -> np.ndarray:
     def get_unitary_and_grad(
         self,
         params: RealVector = [],
-    ) -> tuple[UnitaryMatrix, np.ndarray]:
+    ) -> tuple[UnitaryMatrix, npt.NDArray[np.complex128]]:
         """
         Return the unitary and gradient for this gate.
 
@@ -115,7 +116,7 @@ def get_unitary_and_grad(
         utry, grads = self.gate.get_unitary_and_grad(f_params)  # type: ignore
         return utry, grads[self.unfixed_param_idxs, :, :]
 
-    def optimize(self, env_matrix: np.ndarray) -> list[float]:
+    def optimize(self, env_matrix: npt.NDArray[np.complex128]) -> list[float]:
         """
         Return the optimal parameters with respect to an environment matrix.
 
@@ -127,6 +128,16 @@ def optimize(self, env_matrix: np.ndarray) -> list[float]:
             if i in self.unfixed_param_idxs
         ]
 
+    def __eq__(self, other: object) -> bool:
+        return (
+            isinstance(other, FrozenParameterGate)
+            and self.gate == other.gate
+            and self.frozen_params == other.frozen_params
+        )
+
+    def __hash__(self) -> int:
+        return hash((self.gate, tuple(self.frozen_params.items())))
+
     @property
     def qasm_name(self) -> str:
         """The qasm command for this gate, see Gate for more info."""
diff --git a/bqskit/ir/gates/composed/tagged.py b/bqskit/ir/gates/composed/tagged.py
index 139198cb5..b862e3a8d 100644
--- a/bqskit/ir/gates/composed/tagged.py
+++ b/bqskit/ir/gates/composed/tagged.py
@@ -4,6 +4,7 @@
 from typing import Any
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.ir.gate import Gate
 from bqskit.ir.gates.composedgate import ComposedGate
@@ -55,7 +56,7 @@ def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
 
         return self.gate.get_unitary(params)
 
-    def get_grad(self, params: RealVector = []) -> np.ndarray:
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
         """
         Return the gradient for this gate.
 
@@ -69,7 +70,7 @@ def get_grad(self, params: RealVector = []) -> np.ndarray:
     def get_unitary_and_grad(
         self,
         params: RealVector = [],
-    ) -> tuple[UnitaryMatrix, np.ndarray]:
+    ) -> tuple[UnitaryMatrix, npt.NDArray[np.complex128]]:
         """
         Return the unitary and gradient for this gate.
 
@@ -80,7 +81,7 @@ def get_unitary_and_grad(
 
         return self.gate.get_unitary_and_grad(params)  # type: ignore
 
-    def optimize(self, env_matrix: np.ndarray) -> list[float]:
+    def optimize(self, env_matrix: npt.NDArray[np.complex128]) -> list[float]:
         """
         Return the optimal parameters with respect to an environment matrix.
 
@@ -91,3 +92,13 @@ def optimize(self, env_matrix: np.ndarray) -> list[float]:
 
         self.check_env_matrix(env_matrix)
         return self.gate.optimize(env_matrix)  # type: ignore
+
+    def __eq__(self, other: object) -> bool:
+        return (
+            isinstance(other, TaggedGate)
+            and self.gate == other.gate
+            and self.tag == other.tag
+        )
+
+    def __hash__(self) -> int:
+        return hash((self.gate, self.tag))
diff --git a/bqskit/ir/gates/composed/vlg.py b/bqskit/ir/gates/composed/vlg.py
index b878e92fb..e692e407d 100644
--- a/bqskit/ir/gates/composed/vlg.py
+++ b/bqskit/ir/gates/composed/vlg.py
@@ -5,6 +5,7 @@
 from typing import Sequence
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.ir.gate import Gate
 from bqskit.ir.gates.composedgate import ComposedGate
@@ -74,7 +75,8 @@ def __init__(
             raise ValueError('VLGs require at least 1 locations.')
 
         self.gate = gate
-        self._name = 'VariableLocationGate(%s)' % gate.name
+        name_str_tuple = (gate.name, locations, radixes)
+        self._name = 'VariableLocationGate(%s, %s, %s)' % name_str_tuple
         self.locations = list(locations)
         self._num_qudits = len(set(sum((tuple(l) for l in locations), tuple())))
 
@@ -127,7 +129,7 @@ def get_location(self, params: RealVector) -> tuple[int, ...]:
     def split_params(
         self,
         params: RealVector,
-    ) -> tuple[np.ndarray, np.ndarray]:
+    ) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
         """Split params into subgate params and location params."""
         return (
             np.array(params[:self.gate.num_params]),
@@ -145,7 +147,7 @@ def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
         PGPT = P @ np.kron(G, self.I) @ P.T
         return UnitaryMatrix.closest_to(PGPT, self.radixes)
 
-    def get_grad(self, params: RealVector = []) -> np.ndarray:
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
         """
         Return the gradient for this gate.
 
@@ -156,7 +158,7 @@ def get_grad(self, params: RealVector = []) -> np.ndarray:
     def get_unitary_and_grad(
         self,
         params: RealVector = [],
-    ) -> tuple[UnitaryMatrix, np.ndarray]:
+    ) -> tuple[UnitaryMatrix, npt.NDArray[np.complex128]]:
         """
         Return the unitary and gradient for this gate.
 
@@ -183,10 +185,21 @@ def get_unitary_and_grad(
         U = UnitaryMatrix.closest_to(PGPT, self.radixes)
         return U, np.concatenate([dG, dP])
 
-    def optimize(self, env_matrix: np.ndarray) -> list[float]:
+    def optimize(self, env_matrix: npt.NDArray[np.complex128]) -> list[float]:
         """
         Return the optimal parameters with respect to an environment matrix.
 
         See :class:`LocallyOptimizableUnitary` for more info.
         """
         raise NotImplementedError()
+
+    def __eq__(self, other: object) -> bool:
+        return (
+            isinstance(other, VariableLocationGate)
+            and self.gate == other.gate
+            and self.locations == other.locations
+            and self.radixes == other.radixes
+        )
+
+    def __hash__(self) -> int:
+        return hash((self.gate, tuple(self.locations), self.radixes))
diff --git a/bqskit/ir/gates/constant/__init__.py b/bqskit/ir/gates/constant/__init__.py
index 502af7a76..83c264041 100644
--- a/bqskit/ir/gates/constant/__init__.py
+++ b/bqskit/ir/gates/constant/__init__.py
@@ -15,6 +15,7 @@
 from bqskit.ir.gates.constant.s import SGate
 from bqskit.ir.gates.constant.sdg import SdgGate
 from bqskit.ir.gates.constant.sqrtcnot import SqrtCNOTGate
+from bqskit.ir.gates.constant.sqrtiswap import SqrtISwapGate
 from bqskit.ir.gates.constant.swap import SwapGate
 from bqskit.ir.gates.constant.sx import SqrtXGate
 from bqskit.ir.gates.constant.sx import SXGate
@@ -45,6 +46,7 @@
     'SqrtCNOTGate',
     'SwapGate',
     'SqrtXGate',
+    'SqrtISwapGate',
     'SXGate',
     'TGate',
     'TdgGate',
diff --git a/bqskit/ir/gates/constant/identity.py b/bqskit/ir/gates/constant/identity.py
index 5c0260596..162408fd9 100644
--- a/bqskit/ir/gates/constant/identity.py
+++ b/bqskit/ir/gates/constant/identity.py
@@ -50,3 +50,13 @@ def get_qasm_gate_def(self) -> str:
         body_stmts = ['\tU(0,0,0) %s;' % sym for sym in param_symbols]
         body = '\n'.join(body_stmts)
         return f'{header}\n{{\n{body}\n}}\n'
+
+    def __eq__(self, other: object) -> bool:
+        return (
+            isinstance(other, IdentityGate)
+            and self.num_qudits == other.num_qudits
+            and self.radixes == other.radixes
+        )
+
+    def __hash__(self) -> int:
+        return hash((self.num_qudits, self.radixes))
diff --git a/bqskit/ir/gates/constant/sqrtiswap.py b/bqskit/ir/gates/constant/sqrtiswap.py
new file mode 100644
index 000000000..62e688e0f
--- /dev/null
+++ b/bqskit/ir/gates/constant/sqrtiswap.py
@@ -0,0 +1,35 @@
+"""This module implements the SqrtISwapGate."""
+from __future__ import annotations
+
+import math
+
+from bqskit.ir.gates.constantgate import ConstantGate
+from bqskit.ir.gates.qubitgate import QubitGate
+from bqskit.qis.unitary.unitarymatrix import UnitaryMatrix
+
+
+class SqrtISwapGate(ConstantGate, QubitGate):
+    """
+    The square root two qubit swap and phase iSWAP gate.
+
+    The SqrtISwap gate is given by the following unitary:
+
+    .. math::
+
+        \\begin{pmatrix}
+        1 & 0 & 0 & 0 \\\\
+        0 & \\frac{1}{\\sqrt{2}} & \\frac{i}{\\sqrt{2}} & 0 \\\\
+        0 & \\frac{i}{\\sqrt{2}} & \\frac{1}{\\sqrt{2}} & 0 \\\\
+        0 & 0 & 0 & 1 \\\\
+        \\end{pmatrix}
+    """
+
+    _num_qudits = 2
+    _utry = UnitaryMatrix(
+        [
+            [1, 0, 0, 0],
+            [0, 1 / math.sqrt(2), 1j / math.sqrt(2), 0],
+            [0, 1j / math.sqrt(2), 1 / math.sqrt(2), 0],
+            [0, 0, 0, 1],
+        ],
+    )
diff --git a/bqskit/ir/gates/constant/unitary.py b/bqskit/ir/gates/constant/unitary.py
index 316d332be..fb3f80c50 100644
--- a/bqskit/ir/gates/constant/unitary.py
+++ b/bqskit/ir/gates/constant/unitary.py
@@ -28,3 +28,12 @@ def __init__(
         self._utry = UnitaryMatrix(utry, radixes)
         self._num_qudits = self._utry.num_qudits
         self._radixes = self._utry.radixes
+
+    def __eq__(self, other: object) -> bool:
+        return (
+            isinstance(other, ConstantUnitaryGate)
+            and self._utry == other._utry
+        )
+
+    def __hash__(self) -> int:
+        return hash(self._utry)
diff --git a/bqskit/ir/gates/constantgate.py b/bqskit/ir/gates/constantgate.py
index 59981033a..c6edaccc7 100644
--- a/bqskit/ir/gates/constantgate.py
+++ b/bqskit/ir/gates/constantgate.py
@@ -2,6 +2,7 @@
 from __future__ import annotations
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.ir.gate import Gate
 from bqskit.qis.unitary.differentiable import DifferentiableUnitary
@@ -27,7 +28,7 @@ def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
         self.check_parameters(params)
         return getattr(self, '_utry')
 
-    def get_grad(self, params: RealVector = []) -> np.ndarray:
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
         """
         Return the gradient for this gate.
 
@@ -39,7 +40,7 @@ def get_grad(self, params: RealVector = []) -> np.ndarray:
     def get_unitary_and_grad(
         self,
         params: RealVector = [],
-    ) -> tuple[UnitaryMatrix, np.ndarray]:
+    ) -> tuple[UnitaryMatrix, npt.NDArray[np.complex128]]:
         """
         Return the unitary and gradient for this gate.
 
@@ -48,7 +49,7 @@ def get_unitary_and_grad(
         self.check_parameters(params)
         return getattr(self, '_utry'), np.array([])
 
-    def optimize(self, env_matrix: np.ndarray) -> list[float]:
+    def optimize(self, env_matrix: npt.NDArray[np.complex128]) -> list[float]:
         """
         Return the optimal parameters with respect to an environment matrix.
 
diff --git a/bqskit/ir/gates/parameterized/__init__.py b/bqskit/ir/gates/parameterized/__init__.py
index dd24ee782..aaeb56e88 100644
--- a/bqskit/ir/gates/parameterized/__init__.py
+++ b/bqskit/ir/gates/parameterized/__init__.py
@@ -1,6 +1,9 @@
 """This package contains parameterized gates."""
 from __future__ import annotations
 
+from bqskit.ir.gates.parameterized.crx import CRXGate
+from bqskit.ir.gates.parameterized.cry import CRYGate
+from bqskit.ir.gates.parameterized.crz import CRZGate
 from bqskit.ir.gates.parameterized.pauli import PauliGate
 from bqskit.ir.gates.parameterized.rx import RXGate
 from bqskit.ir.gates.parameterized.rxx import RXXGate
@@ -15,6 +18,9 @@
 from bqskit.ir.gates.parameterized.unitary import VariableUnitaryGate
 
 __all__ = [
+    'CRXGate',
+    'CRYGate',
+    'CRZGate',
     'PauliGate',
     'RXGate',
     'RXXGate',
diff --git a/bqskit/ir/gates/parameterized/crx.py b/bqskit/ir/gates/parameterized/crx.py
new file mode 100644
index 000000000..0f49fcbdc
--- /dev/null
+++ b/bqskit/ir/gates/parameterized/crx.py
@@ -0,0 +1,74 @@
+"""This module implements the CRXGate."""
+from __future__ import annotations
+
+import numpy as np
+import numpy.typing as npt
+
+from bqskit.ir.gates.qubitgate import QubitGate
+from bqskit.qis.unitary.differentiable import DifferentiableUnitary
+from bqskit.qis.unitary.unitary import RealVector
+from bqskit.qis.unitary.unitarymatrix import UnitaryMatrix
+from bqskit.utils.cachedclass import CachedClass
+
+
+class CRXGate(
+    QubitGate,
+    DifferentiableUnitary,
+    CachedClass,
+):
+    """
+    A gate representing a controlled X rotation.
+
+    It is given by the following parameterized unitary:
+
+    .. math::
+
+        \\begin{pmatrix}
+        1 & 0 & 0 & 0 \\\\
+        0 & 1 & 0 & 0 \\\\
+        0 & 0 & \\cos{\\frac{\\theta}{2}} & -\\sin{\\frac{\\theta}{2}}i \\\\
+        0 & 0 & -\\sin{\\frac{\\theta}{2}}i & \\cos{\\frac{\\theta}{2}} \\\\
+        \\end{pmatrix}
+    """
+
+    _num_qudits = 2
+    _num_params = 1
+    _qasm_name = 'crx'
+
+    def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
+        """Return the unitary for this gate, see :class:`Unitary` for more."""
+        self.check_parameters(params)
+
+        cos = np.cos(params[0] / 2)
+        sin = -1j * np.sin(params[0] / 2)
+
+        return UnitaryMatrix(
+            [
+                [1, 0, 0, 0],
+                [0, 1, 0, 0],
+                [0, 0, cos, sin],
+                [0, 0, sin, cos],
+            ],
+        )
+
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
+        """
+        Return the gradient for this gate.
+
+        See :class:`DifferentiableUnitary` for more info.
+        """
+        self.check_parameters(params)
+
+        dcos = -np.sin(params[0] / 2) / 2
+        dsin = -1j * np.cos(params[0] / 2) / 2
+
+        return np.array(
+            [
+                [
+                    [0, 0, 0, 0],
+                    [0, 0, 0, 0],
+                    [0, 0, dcos, dsin],
+                    [0, 0, dsin, dcos],
+                ],
+            ], dtype=np.complex128,
+        )
diff --git a/bqskit/ir/gates/parameterized/cry.py b/bqskit/ir/gates/parameterized/cry.py
new file mode 100644
index 000000000..871735fc7
--- /dev/null
+++ b/bqskit/ir/gates/parameterized/cry.py
@@ -0,0 +1,74 @@
+"""This module implements the CRYGate."""
+from __future__ import annotations
+
+import numpy as np
+import numpy.typing as npt
+
+from bqskit.ir.gates.qubitgate import QubitGate
+from bqskit.qis.unitary.differentiable import DifferentiableUnitary
+from bqskit.qis.unitary.unitary import RealVector
+from bqskit.qis.unitary.unitarymatrix import UnitaryMatrix
+from bqskit.utils.cachedclass import CachedClass
+
+
+class CRYGate(
+    QubitGate,
+    DifferentiableUnitary,
+    CachedClass,
+):
+    """
+    A gate representing a controlled Y rotation.
+
+    It is given by the following parameterized unitary:
+
+    .. math::
+
+        \\begin{pmatrix}
+        1 & 0 & 0 & 0 \\\\
+        0 & 1 & 0 & 0 \\\\
+        0 & 0 & \\cos{\\frac{\\theta}{2}} & -\\sin{\\frac{\\theta}{2}} \\\\
+        0 & 0 & \\sin{\\frac{\\theta}{2}} & \\cos{\\frac{\\theta}{2}} \\\\
+        \\end{pmatrix}
+    """
+
+    _num_qudits = 2
+    _num_params = 1
+    _qasm_name = 'cry'
+
+    def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
+        """Return the unitary for this gate, see :class:`Unitary` for more."""
+        self.check_parameters(params)
+
+        cos = np.cos(params[0] / 2)
+        sin = np.sin(params[0] / 2)
+
+        return UnitaryMatrix(
+            [
+                [1, 0, 0, 0],
+                [0, 1, 0, 0],
+                [0, 0, cos, -sin],
+                [0, 0, sin, cos],
+            ],
+        )
+
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
+        """
+        Return the gradient for this gate.
+
+        See :class:`DifferentiableUnitary` for more info.
+        """
+        self.check_parameters(params)
+
+        dcos = -np.sin(params[0] / 2) / 2
+        dsin = -1j * np.cos(params[0] / 2) / 2
+
+        return np.array(
+            [
+                [
+                    [0, 0, 0, 0],
+                    [0, 0, 0, 0],
+                    [0, 0, dcos, -dsin],
+                    [0, 0, dsin, dcos],
+                ],
+            ], dtype=np.complex128,
+        )
diff --git a/bqskit/ir/gates/parameterized/crz.py b/bqskit/ir/gates/parameterized/crz.py
new file mode 100644
index 000000000..0fc5b6a21
--- /dev/null
+++ b/bqskit/ir/gates/parameterized/crz.py
@@ -0,0 +1,74 @@
+"""This module implements the CRZGate."""
+from __future__ import annotations
+
+import numpy as np
+import numpy.typing as npt
+
+from bqskit.ir.gates.qubitgate import QubitGate
+from bqskit.qis.unitary.differentiable import DifferentiableUnitary
+from bqskit.qis.unitary.unitary import RealVector
+from bqskit.qis.unitary.unitarymatrix import UnitaryMatrix
+from bqskit.utils.cachedclass import CachedClass
+
+
+class CRZGate(
+    QubitGate,
+    DifferentiableUnitary,
+    CachedClass,
+):
+    """
+    A gate representing a controlled Z rotation.
+
+    It is given by the following parameterized unitary:
+
+    .. math::
+
+        \\begin{pmatrix}
+        1 & 0 & 0 & 0 \\\\
+        0 & 1 & 0 & 0 \\\\
+        0 & 0 & \\exp({-i\\frac{\\theta}{2}}) & 0 \\\\
+        0 & 0 & 0 & \\exp({i\\frac{\\theta}{2}}) \\\\
+        \\end{pmatrix}
+    """
+
+    _num_qudits = 2
+    _num_params = 1
+    _qasm_name = 'crz'
+
+    def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
+        """Return the unitary for this gate, see :class:`Unitary` for more."""
+        self.check_parameters(params)
+
+        pos = np.exp(1j * params[0] / 2)
+        neg = np.exp(-1j * params[0] / 2)
+
+        return UnitaryMatrix(
+            [
+                [1, 0, 0, 0],
+                [0, 1, 0, 0],
+                [0, 0, neg, 0],
+                [0, 0, 0, pos],
+            ],
+        )
+
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
+        """
+        Return the gradient for this gate.
+
+        See :class:`DifferentiableUnitary` for more info.
+        """
+        self.check_parameters(params)
+
+        dpos = 1j / 2 * np.exp(1j * params[0] / 2)
+        dneg = -1j / 2 * np.exp(-1j * params[0] / 2)
+
+        return np.array(
+            [
+                [
+                    [0, 0, 0, 0],
+                    [0, 0, 0, 0],
+                    [0, 0, dneg, 0],
+                    [0, 0, 0, dpos],
+                ],
+            ], dtype=np.complex128,
+        )
diff --git a/bqskit/ir/gates/parameterized/pauli.py b/bqskit/ir/gates/parameterized/pauli.py
index 482d8aac8..d6d900bb5 100644
--- a/bqskit/ir/gates/parameterized/pauli.py
+++ b/bqskit/ir/gates/parameterized/pauli.py
@@ -1,9 +1,10 @@
 """This module implements the PauliGate."""
 from __future__ import annotations
 
-import os
+from typing import Any
 
 import numpy as np
+import numpy.typing as npt
 import scipy as sp
 
 from bqskit.ir.gates.qubitgate import QubitGate
@@ -12,6 +13,7 @@
 from bqskit.qis.unitary.optimizable import LocallyOptimizableUnitary
 from bqskit.qis.unitary.unitary import RealVector
 from bqskit.qis.unitary.unitarymatrix import UnitaryMatrix
+from bqskit.utils.docs import building_docs
 from bqskit.utils.math import dexpmv
 from bqskit.utils.math import dot_product
 from bqskit.utils.math import pauli_expansion
@@ -47,11 +49,12 @@ def __init__(self, num_qudits: int) -> None:
         if num_qudits <= 0:
             raise ValueError('Expected positive integer, got %d' % num_qudits)
 
+        self._name = 'PauliGate(%s)' % num_qudits
         self._num_qudits = num_qudits
         self.paulis = PauliMatrices(self.num_qudits)
         self._num_params = len(self.paulis)
-        if 'READTHEDOCS' in os.environ:
-            self.sigmav = np.array([])
+        if building_docs():
+            self.sigmav: npt.NDArray[Any] = np.array([])
         else:
             self.sigmav = (-1j / 2) * self.paulis.numpy
 
@@ -63,7 +66,7 @@ def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
         eiH = sp.linalg.expm(H)
         return UnitaryMatrix(eiH, check_arguments=False)
 
-    def get_grad(self, params: RealVector = []) -> np.ndarray:
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
         """
         Return the gradient for this gate.
 
@@ -78,7 +81,7 @@ def get_grad(self, params: RealVector = []) -> np.ndarray:
     def get_unitary_and_grad(
         self,
         params: RealVector = [],
-    ) -> tuple[UnitaryMatrix, np.ndarray]:
+    ) -> tuple[UnitaryMatrix, npt.NDArray[np.complex128]]:
         """
         Return the unitary and gradient for this gate.
 
@@ -90,7 +93,7 @@ def get_unitary_and_grad(
         U, dU = dexpmv(H, self.sigmav)
         return UnitaryMatrix(U, check_arguments=False), dU
 
-    def optimize(self, env_matrix: np.ndarray) -> list[float]:
+    def optimize(self, env_matrix: npt.NDArray[np.complex128]) -> list[float]:
         """
         Return the optimal parameters with respect to an environment matrix.
 
diff --git a/bqskit/ir/gates/parameterized/rx.py b/bqskit/ir/gates/parameterized/rx.py
index 640c4e291..7db98bf0b 100644
--- a/bqskit/ir/gates/parameterized/rx.py
+++ b/bqskit/ir/gates/parameterized/rx.py
@@ -2,6 +2,7 @@
 from __future__ import annotations
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.ir.gates.qubitgate import QubitGate
 from bqskit.qis.unitary.differentiable import DifferentiableUnitary
@@ -48,7 +49,7 @@ def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
             ],
         )
 
-    def get_grad(self, params: RealVector = []) -> np.ndarray:
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
         """
         Return the gradient for this gate.
 
@@ -68,7 +69,7 @@ def get_grad(self, params: RealVector = []) -> np.ndarray:
             ], dtype=np.complex128,
         )
 
-    def optimize(self, env_matrix: np.ndarray) -> list[float]:
+    def optimize(self, env_matrix: npt.NDArray[np.complex128]) -> list[float]:
         """
         Return the optimal parameters with respect to an environment matrix.
 
diff --git a/bqskit/ir/gates/parameterized/rxx.py b/bqskit/ir/gates/parameterized/rxx.py
index 16dd40b61..7696ed0e3 100644
--- a/bqskit/ir/gates/parameterized/rxx.py
+++ b/bqskit/ir/gates/parameterized/rxx.py
@@ -2,6 +2,7 @@
 from __future__ import annotations
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.ir.gates.qubitgate import QubitGate
 from bqskit.qis.unitary.differentiable import DifferentiableUnitary
@@ -52,7 +53,7 @@ def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
             ],
         )
 
-    def get_grad(self, params: RealVector = []) -> np.ndarray:
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
         """
         Return the gradient for this gate.
 
@@ -74,7 +75,7 @@ def get_grad(self, params: RealVector = []) -> np.ndarray:
             ], dtype=np.complex128,
         )
 
-    def optimize(self, env_matrix: np.ndarray) -> list[float]:
+    def optimize(self, env_matrix: npt.NDArray[np.complex128]) -> list[float]:
         """
         Return the optimal parameters with respect to an environment matrix.
 
diff --git a/bqskit/ir/gates/parameterized/ry.py b/bqskit/ir/gates/parameterized/ry.py
index 6f53071e8..c094aa355 100644
--- a/bqskit/ir/gates/parameterized/ry.py
+++ b/bqskit/ir/gates/parameterized/ry.py
@@ -2,6 +2,7 @@
 from __future__ import annotations
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.ir.gates.qubitgate import QubitGate
 from bqskit.qis.unitary.differentiable import DifferentiableUnitary
@@ -48,7 +49,7 @@ def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
             ],
         )
 
-    def get_grad(self, params: RealVector = []) -> np.ndarray:
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
         """
         Return the gradient for this gate.
 
@@ -68,7 +69,7 @@ def get_grad(self, params: RealVector = []) -> np.ndarray:
             ], dtype=np.complex128,
         )
 
-    def optimize(self, env_matrix: np.ndarray) -> list[float]:
+    def optimize(self, env_matrix: npt.NDArray[np.complex128]) -> list[float]:
         """
         Return the optimal parameters with respect to an environment matrix.
 
diff --git a/bqskit/ir/gates/parameterized/ryy.py b/bqskit/ir/gates/parameterized/ryy.py
index 7e5512eb4..dd9e7440d 100644
--- a/bqskit/ir/gates/parameterized/ryy.py
+++ b/bqskit/ir/gates/parameterized/ryy.py
@@ -2,6 +2,7 @@
 from __future__ import annotations
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.ir.gates.qubitgate import QubitGate
 from bqskit.qis.unitary.differentiable import DifferentiableUnitary
@@ -51,7 +52,7 @@ def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
             ],
         )
 
-    def get_grad(self, params: RealVector = []) -> np.ndarray:
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
         """
         Return the gradient for this gate.
 
diff --git a/bqskit/ir/gates/parameterized/rz.py b/bqskit/ir/gates/parameterized/rz.py
index 0df8372df..c7a0eba5e 100644
--- a/bqskit/ir/gates/parameterized/rz.py
+++ b/bqskit/ir/gates/parameterized/rz.py
@@ -2,6 +2,7 @@
 from __future__ import annotations
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.ir.gates.qubitgate import QubitGate
 from bqskit.qis.unitary.differentiable import DifferentiableUnitary
@@ -42,7 +43,7 @@ def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
             ],
         )
 
-    def get_grad(self, params: RealVector = []) -> np.ndarray:
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
         """
         Return the gradient for this gate.
 
diff --git a/bqskit/ir/gates/parameterized/rzz.py b/bqskit/ir/gates/parameterized/rzz.py
index 892f7696d..967bef2d0 100644
--- a/bqskit/ir/gates/parameterized/rzz.py
+++ b/bqskit/ir/gates/parameterized/rzz.py
@@ -2,6 +2,7 @@
 from __future__ import annotations
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.ir.gates.qubitgate import QubitGate
 from bqskit.qis.unitary.differentiable import DifferentiableUnitary
@@ -50,7 +51,7 @@ def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
             ],
         )
 
-    def get_grad(self, params: RealVector = []) -> np.ndarray:
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
         """
         Return the gradient for this gate.
 
diff --git a/bqskit/ir/gates/parameterized/u1.py b/bqskit/ir/gates/parameterized/u1.py
index f28b02b52..57a0cbade 100644
--- a/bqskit/ir/gates/parameterized/u1.py
+++ b/bqskit/ir/gates/parameterized/u1.py
@@ -2,6 +2,7 @@
 from __future__ import annotations
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.ir.gates.qubitgate import QubitGate
 from bqskit.qis.unitary.differentiable import DifferentiableUnitary
@@ -47,7 +48,7 @@ def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
             ],
         )
 
-    def get_grad(self, params: RealVector = []) -> np.ndarray:
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
         """
         Return the gradient for this gate.
 
@@ -66,7 +67,7 @@ def get_grad(self, params: RealVector = []) -> np.ndarray:
             ], dtype=np.complex128,
         )
 
-    def optimize(self, env_matrix: np.ndarray) -> list[float]:
+    def optimize(self, env_matrix: npt.NDArray[np.complex128]) -> list[float]:
         """
         Return the optimal parameters with respect to an environment matrix.
 
diff --git a/bqskit/ir/gates/parameterized/u2.py b/bqskit/ir/gates/parameterized/u2.py
index 11925941d..9f6a3816e 100644
--- a/bqskit/ir/gates/parameterized/u2.py
+++ b/bqskit/ir/gates/parameterized/u2.py
@@ -2,6 +2,7 @@
 from __future__ import annotations
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.ir.gates.qubitgate import QubitGate
 from bqskit.qis.unitary.differentiable import DifferentiableUnitary
@@ -44,7 +45,7 @@ def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
             ],
         )
 
-    def get_grad(self, params: RealVector = []) -> np.ndarray:
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
         """
         Return the gradient for this gate.
 
diff --git a/bqskit/ir/gates/parameterized/u3.py b/bqskit/ir/gates/parameterized/u3.py
index 311f79cad..97d2fb964 100644
--- a/bqskit/ir/gates/parameterized/u3.py
+++ b/bqskit/ir/gates/parameterized/u3.py
@@ -2,6 +2,7 @@
 from __future__ import annotations
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.ir.gates.qubitgate import QubitGate
 from bqskit.qis.unitary.differentiable import DifferentiableUnitary
@@ -50,7 +51,7 @@ def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
             ],
         )
 
-    def get_grad(self, params: RealVector = []) -> np.ndarray:
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
         """
         Return the gradient for this gate.
 
diff --git a/bqskit/ir/gates/parameterized/u8.py b/bqskit/ir/gates/parameterized/u8.py
index a53a4cecd..af08b1a86 100644
--- a/bqskit/ir/gates/parameterized/u8.py
+++ b/bqskit/ir/gates/parameterized/u8.py
@@ -2,6 +2,7 @@
 from __future__ import annotations
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.ir.gates.qutritgate import QutritGate
 from bqskit.qis.unitary.differentiable import DifferentiableUnitary
@@ -58,7 +59,7 @@ def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
             ],
         )
 
-    def get_grad(self, params: RealVector = []) -> np.ndarray:
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
         """
         Return the gradient for this gate.
 
diff --git a/bqskit/ir/gates/parameterized/unitary.py b/bqskit/ir/gates/parameterized/unitary.py
index 39af713dc..e571df0e4 100644
--- a/bqskit/ir/gates/parameterized/unitary.py
+++ b/bqskit/ir/gates/parameterized/unitary.py
@@ -4,11 +4,13 @@
 from typing import Sequence
 
 import numpy as np
+import numpy.typing as npt
 import scipy as sp
 
 from bqskit.ir.gate import Gate
 from bqskit.qis.unitary.optimizable import LocallyOptimizableUnitary
 from bqskit.qis.unitary.unitary import RealVector
+from bqskit.qis.unitary.unitarymatrix import UnitaryLike
 from bqskit.qis.unitary.unitarymatrix import UnitaryMatrix
 from bqskit.utils.typing import is_valid_radixes
 
@@ -38,12 +40,14 @@ def __init__(self, num_qudits: int, radixes: Sequence[int] = []) -> None:
         if not is_valid_radixes(radixes, num_qudits):
             raise TypeError('Invalid radixes.')
 
-        self._num_qudits = num_qudits
+        self._num_qudits = int(num_qudits)
         self._radixes = tuple(radixes)
         self._dim = int(np.prod(self.radixes))
         self.shape = (self.dim, self.dim)
         self._num_params = 2 * self.dim**2
-        self._name = 'VariableUnitaryGate(%d)' % self.num_qudits
+        self._name = 'VariableUnitaryGate(%d, %s)' % (
+            self.num_qudits, str(self.radixes),
+        )
 
     def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
         """
@@ -61,7 +65,7 @@ def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
         x = real + imag
         return UnitaryMatrix.closest_to(np.reshape(x, self.shape), self.radixes)
 
-    def optimize(self, env_matrix: np.ndarray) -> list[float]:
+    def optimize(self, env_matrix: npt.NDArray[np.complex128]) -> list[float]:
         """
         Return the optimal parameters with respect to an environment matrix.
 
@@ -71,3 +75,21 @@ def optimize(self, env_matrix: np.ndarray) -> list[float]:
         U, _, Vh = sp.linalg.svd(env_matrix)
         x = np.reshape(Vh.conj().T @ U.conj().T, (self.num_params // 2,))
         return list(np.real(x)) + list(np.imag(x))
+
+    @staticmethod
+    def get_params(utry: UnitaryLike) -> RealVector:
+        """Return the params for this gate, given a unitary matrix."""
+        num_elems = len(utry) ** 2
+        real = np.reshape(np.real(utry), num_elems)
+        imag = np.reshape(np.imag(utry), num_elems)
+        return np.concatenate([real, imag])
+
+    def __eq__(self, other: object) -> bool:
+        return (
+            isinstance(other, VariableUnitaryGate)
+            and self.num_qudits == other.num_qudits
+            and self.radixes == other.radixes
+        )
+
+    def __hash__(self) -> int:
+        return hash((self.num_qudits, self.radixes))
diff --git a/bqskit/ir/interval.py b/bqskit/ir/interval.py
index 69bc8df1b..53b061fd9 100644
--- a/bqskit/ir/interval.py
+++ b/bqskit/ir/interval.py
@@ -108,10 +108,7 @@ def indices(self) -> list[int]:
 
     def __contains__(self, cycle_index: object) -> bool:
         """Return true if `cycle_index` is inside this interval."""
-        if not is_integer(cycle_index):
-            return False
-
-        return self.lower <= cycle_index <= self.upper
+        return self.lower <= cycle_index <= self.upper  # type: ignore
 
     def __iter__(self) -> Iterator[int]:
         """Return an iterator for all indices contained in the interval."""
diff --git a/bqskit/ir/iterator.py b/bqskit/ir/iterator.py
index eb3e31e28..6e6659215 100644
--- a/bqskit/ir/iterator.py
+++ b/bqskit/ir/iterator.py
@@ -195,13 +195,6 @@ def decrement_iter(self) -> None:
 
     def step(self) -> None:
         """Move the iterator one step."""
-        if (
-            self.num_ops != self.circuit.num_operations
-            or self.num_cycles != self.circuit.num_cycles
-            or self.num_qudits != self.circuit.num_qudits
-        ):
-            raise RuntimeError('Circuit changed under iteration.')
-
         if not self.reverse:
             self.increment_iter()
         else:
@@ -212,6 +205,13 @@ def step(self) -> None:
             raise StopIteration
 
     def __next__(self) -> Operation | tuple[int, Operation]:
+        if (
+            self.num_ops != self.circuit.num_operations
+            or self.num_cycles != self.circuit.num_cycles
+            or self.num_qudits != self.circuit.num_qudits
+        ):
+            raise RuntimeError('Circuit changed under iteration.')
+
         while True:
             self.step()
             op = self.circuit._circuit[self.cycle][self.qudit]
diff --git a/bqskit/ir/location.py b/bqskit/ir/location.py
index 83d142a2e..e114d0573 100644
--- a/bqskit/ir/location.py
+++ b/bqskit/ir/location.py
@@ -3,6 +3,7 @@
 
 import logging
 from typing import Any
+from typing import Iterator
 from typing import overload
 from typing import Sequence
 from typing import Union
@@ -186,6 +187,9 @@ def __repr__(self) -> str:
     def __hash__(self) -> int:
         return self._location.__hash__()
 
+    def __iter__(self) -> Iterator[int]:
+        return self._location.__iter__()
+
     def __eq__(self, other: object) -> bool:
         if isinstance(other, CircuitLocation):
             return self._location.__eq__(other._location)
diff --git a/bqskit/ir/operation.py b/bqskit/ir/operation.py
index fd6bc5aa6..0d9ced7d5 100644
--- a/bqskit/ir/operation.py
+++ b/bqskit/ir/operation.py
@@ -4,6 +4,7 @@
 from typing import Any
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.ir.gate import Gate
 from bqskit.ir.gates import FrozenParameterGate
@@ -109,7 +110,7 @@ def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
 
         return self.gate.get_unitary(self.params)
 
-    def get_grad(self, params: RealVector = []) -> np.ndarray:
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
         """
         Return the gradient for this operation.
 
@@ -123,7 +124,7 @@ def get_grad(self, params: RealVector = []) -> np.ndarray:
     def get_unitary_and_grad(
         self,
         params: RealVector = [],
-    ) -> tuple[UnitaryMatrix, np.ndarray]:
+    ) -> tuple[UnitaryMatrix, npt.NDArray[np.complex128]]:
         """
         Return the unitary and gradient for this gate.
 
diff --git a/bqskit/ir/opt/cost/differentiable.py b/bqskit/ir/opt/cost/differentiable.py
index a395a3f27..aeacbd021 100644
--- a/bqskit/ir/opt/cost/differentiable.py
+++ b/bqskit/ir/opt/cost/differentiable.py
@@ -4,6 +4,7 @@
 import abc
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.ir.opt.cost.function import CostFunction
 from bqskit.ir.opt.cost.residual import ResidualsFunction
@@ -24,13 +25,13 @@ class DifferentiableCostFunction(CostFunction):
     """
 
     @abc.abstractmethod
-    def get_grad(self, params: RealVector) -> np.ndarray:
+    def get_grad(self, params: RealVector) -> npt.NDArray[np.float64]:
         """Return the cost gradient given the input parameters."""
 
     def get_cost_and_grad(
         self,
         params: RealVector,
-    ) -> tuple[float, np.ndarray]:
+    ) -> tuple[float, npt.NDArray[np.float64]]:
         """Return the cost and gradient given the input parameters."""
         return self.get_cost(params), self.get_grad(params)
 
@@ -50,12 +51,12 @@ class DifferentiableResidualsFunction(ResidualsFunction):
     """
 
     @abc.abstractmethod
-    def get_grad(self, params: RealVector) -> np.ndarray:
+    def get_grad(self, params: RealVector) -> npt.NDArray[np.float64]:
         """Return the residuals gradient given the input parameters."""
 
     def get_cost_and_grad(
         self,
         params: RealVector,
-    ) -> tuple[np.ndarray, np.ndarray]:
+    ) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]:
         """Return the residuals and gradient given the input parameters."""
         return self.get_residuals(params), self.get_grad(params)
diff --git a/bqskit/ir/opt/cost/residual.py b/bqskit/ir/opt/cost/residual.py
index 7642a473d..8c19e7718 100644
--- a/bqskit/ir/opt/cost/residual.py
+++ b/bqskit/ir/opt/cost/residual.py
@@ -4,6 +4,7 @@
 import abc
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.ir.opt.cost.function import CostFunction
 from bqskit.qis.unitary.unitary import RealVector
@@ -21,13 +22,13 @@ class ResidualsFunction(CostFunction):
     """
 
     @abc.abstractmethod
-    def get_residuals(self, params: RealVector) -> np.ndarray:
+    def get_residuals(self, params: RealVector) -> npt.NDArray[np.float64]:
         """Return the vector of residuals given the input parameters."""
 
     def __call__(  # type: ignore
         self,
         params: RealVector,
-    ) -> np.ndarray:
+    ) -> npt.NDArray[np.float64]:
         """Return the vector of residuals given the input parameters."""
 
         if not is_sequence(params):
diff --git a/bqskit/ir/opt/instantiater.py b/bqskit/ir/opt/instantiater.py
index c43d007fd..bfe3771b1 100644
--- a/bqskit/ir/opt/instantiater.py
+++ b/bqskit/ir/opt/instantiater.py
@@ -5,6 +5,7 @@
 from typing import TYPE_CHECKING
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.qis.state.state import StateVector
 from bqskit.qis.unitary.unitarymatrix import UnitaryMatrix
@@ -29,8 +30,8 @@ def instantiate(
         self,
         circuit: Circuit,
         target: UnitaryMatrix | StateVector,
-        x0: np.ndarray,
-    ) -> np.ndarray:
+        x0: npt.NDArray[np.float64],
+    ) -> npt.NDArray[np.float64]:
         """
         Instantiate `circuit` to best implement `target`.
 
@@ -76,7 +77,7 @@ def gen_starting_points(
         multistarts: int,
         circuit: Circuit,
         target: UnitaryMatrix | StateVector,
-    ) -> list[np.ndarray]:
+    ) -> list[npt.NDArray[np.float64]]:
         """
         Generate `multistarts` starting points for instantiation.
 
diff --git a/bqskit/ir/opt/instantiaters/minimization.py b/bqskit/ir/opt/instantiaters/minimization.py
index a80c896fd..bab584f34 100644
--- a/bqskit/ir/opt/instantiaters/minimization.py
+++ b/bqskit/ir/opt/instantiaters/minimization.py
@@ -5,6 +5,7 @@
 from typing import TYPE_CHECKING
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.ir.opt.cost.functions import HilbertSchmidtResidualsGenerator
 from bqskit.ir.opt.cost.generator import CostFunctionGenerator
@@ -58,8 +59,8 @@ def instantiate(
         self,
         circuit: Circuit,
         target: UnitaryMatrix | StateVector,
-        x0: np.ndarray,
-    ) -> np.ndarray:
+        x0: npt.NDArray[np.float64],
+    ) -> npt.NDArray[np.float64]:
         """Instantiate `circuit`, see Instantiater for more info."""
         cost = self.cost_fn_gen.gen_cost(circuit, target)
         return self.minimizer.minimize(cost, x0)
diff --git a/bqskit/ir/opt/instantiaters/qfactor.py b/bqskit/ir/opt/instantiaters/qfactor.py
index e5900efed..44e40a7c2 100644
--- a/bqskit/ir/opt/instantiaters/qfactor.py
+++ b/bqskit/ir/opt/instantiaters/qfactor.py
@@ -6,6 +6,7 @@
 from typing import TYPE_CHECKING
 
 import numpy as np
+import numpy.typing as npt
 from bqskitrs import QFactorInstantiatorNative
 
 from bqskit.ir.opt.instantiater import Instantiater
@@ -32,8 +33,8 @@ def instantiate(
         self,
         circuit: Circuit,
         target: UnitaryMatrix | StateVector,
-        x0: np.ndarray,
-    ) -> np.ndarray:
+        x0: npt.NDArray[np.float64],
+    ) -> npt.NDArray[np.float64]:
         """Instantiate `circuit`, see Instantiater for more info."""
         return super().instantiate(circuit, target, x0)
 
diff --git a/bqskit/ir/opt/minimizer.py b/bqskit/ir/opt/minimizer.py
index f9ee52831..c19a7cadc 100644
--- a/bqskit/ir/opt/minimizer.py
+++ b/bqskit/ir/opt/minimizer.py
@@ -5,6 +5,9 @@
 from typing import TYPE_CHECKING
 
 import numpy as np
+import numpy.typing as npt
+
+from bqskit.qis.unitary.unitary import RealVector
 
 if TYPE_CHECKING:
     from bqskit.ir.opt.cost.function import CostFunction
@@ -19,7 +22,10 @@ class Minimizer(abc.ABC):
     """
 
     @abc.abstractmethod
-    def minimize(self, cost: CostFunction, x0: np.ndarray) -> np.ndarray:
+    def minimize(
+        self, cost: CostFunction,
+        x0: RealVector,
+    ) -> npt.NDArray[np.float64]:
         """
         Minimize `cost` starting from the initial point `x0`.
 
diff --git a/bqskit/passes/__init__.py b/bqskit/passes/__init__.py
index 659bc3e39..7ec2e570a 100644
--- a/bqskit/passes/__init__.py
+++ b/bqskit/passes/__init__.py
@@ -28,7 +28,6 @@
     :recursive:
 
     LEAPSynthesisPass
-    OptimizedLEAPPass
     QSearchSynthesisPass
     QFASTDecompositionPass
     QPredictDecompositionPass
@@ -41,6 +40,7 @@
 
     ScanningGateRemovalPass
     WindowOptimizationPass
+    SubstitutePass
 
 .. rubric:: Control Passes
 
@@ -55,6 +55,8 @@
     ForEachBlockPass
     IfThenElsePass
     WhileLoopPass
+    DoThenDecide
+    ParallelDo
 
 .. rubric:: Predicates
 
@@ -89,6 +91,7 @@
     UnfoldPass
     UpdateDataPass
     ToU3Pass
+    BlockConversionPass
 
 .. rubric:: Search Heuristics
 
@@ -110,12 +113,15 @@
     LayerGenerator
     SimpleLayerGenerator
     SeedLayerGenerator
+    StairLayerGenerator
 """
 from __future__ import annotations
 
+from bqskit.passes.control.dothendecide import DoThenDecide
 from bqskit.passes.control.dowhileloop import DoWhileLoopPass
 from bqskit.passes.control.foreach import ForEachBlockPass
 from bqskit.passes.control.ifthenelse import IfThenElsePass
+from bqskit.passes.control.paralleldo import ParallelDo
 from bqskit.passes.control.predicate import PassPredicate
 from bqskit.passes.control.predicates.change import ChangePredicate
 from bqskit.passes.control.predicates.count import GateCountPredicate
@@ -127,22 +133,24 @@
 from bqskit.passes.partitioning.quick import QuickPartitioner
 from bqskit.passes.partitioning.scan import ScanPartitioner
 from bqskit.passes.processing.scan import ScanningGateRemovalPass
+from bqskit.passes.processing.substitute import SubstitutePass
 from bqskit.passes.processing.window import WindowOptimizationPass
 from bqskit.passes.search.frontier import Frontier
 from bqskit.passes.search.generator import LayerGenerator
 from bqskit.passes.search.generators.seed import SeedLayerGenerator
 from bqskit.passes.search.generators.simple import SimpleLayerGenerator
+from bqskit.passes.search.generators.stair import StairLayerGenerator
 from bqskit.passes.search.heuristic import HeuristicFunction
 from bqskit.passes.search.heuristics.astar import AStarHeuristic
 from bqskit.passes.search.heuristics.dijkstra import DijkstraHeuristic
 from bqskit.passes.search.heuristics.greedy import GreedyHeuristic
 from bqskit.passes.synthesis.leap import LEAPSynthesisPass
-from bqskit.passes.synthesis.leap import OptimizedLEAPPass
 from bqskit.passes.synthesis.qfast import QFASTDecompositionPass
 from bqskit.passes.synthesis.qpredict import QPredictDecompositionPass
 from bqskit.passes.synthesis.qsearch import QSearchSynthesisPass
 from bqskit.passes.synthesis.synthesis import SynthesisPass
 from bqskit.passes.util.compress import CompressPass
+from bqskit.passes.util.conversion import BlockConversionPass
 from bqskit.passes.util.converttou3 import ToU3Pass
 from bqskit.passes.util.random import SetRandomSeedPass
 from bqskit.passes.util.record import RecordStatsPass
@@ -164,7 +172,6 @@
     'QuickPartitioner',
     'SynthesisPass',
     'LEAPSynthesisPass',
-    'OptimizedLEAPPass',
     'QSearchSynthesisPass',
     'QFASTDecompositionPass',
     'QPredictDecompositionPass',
@@ -185,4 +192,9 @@
     'LayerGenerator',
     'HeuristicFunction',
     'SeedLayerGenerator',
+    'BlockConversionPass',
+    'StairLayerGenerator',
+    'DoThenDecide',
+    'SubstitutePass',
+    'ParallelDo',
 ]
diff --git a/bqskit/passes/control/__init__.py b/bqskit/passes/control/__init__.py
index ed67d1d55..0f96c933b 100644
--- a/bqskit/passes/control/__init__.py
+++ b/bqskit/passes/control/__init__.py
@@ -1,9 +1,11 @@
 """This package defines passes and objects that control pass execution flow."""
 from __future__ import annotations
 
+from bqskit.passes.control.dothendecide import DoThenDecide
 from bqskit.passes.control.dowhileloop import DoWhileLoopPass
 from bqskit.passes.control.foreach import ForEachBlockPass
 from bqskit.passes.control.ifthenelse import IfThenElsePass
+from bqskit.passes.control.paralleldo import ParallelDo
 from bqskit.passes.control.predicate import PassPredicate
 from bqskit.passes.control.predicates.change import ChangePredicate
 from bqskit.passes.control.predicates.count import GateCountPredicate
@@ -19,4 +21,6 @@
     'GateCountPredicate',
     'NotPredicate',
     'WhileLoopPass',
+    'DoThenDecide',
+    'ParallelDo',
 ]
diff --git a/bqskit/passes/control/dothendecide.py b/bqskit/passes/control/dothendecide.py
new file mode 100644
index 000000000..b29a3dbbb
--- /dev/null
+++ b/bqskit/passes/control/dothendecide.py
@@ -0,0 +1,78 @@
+"""This module implements the DoThenDecide class."""
+from __future__ import annotations
+
+import logging
+from typing import Any
+from typing import Callable
+from typing import Sequence
+
+from bqskit.compiler.basepass import BasePass
+from bqskit.ir.circuit import Circuit
+from bqskit.utils.typing import is_sequence
+
+_logger = logging.getLogger(__name__)
+
+
+class DoThenDecide(BasePass):
+    """
+    The DoThenDecide class.
+
+    This is a control pass that executes a pass or sequence of passes and then
+    conditionally accepts it or reverts to the original state.
+    """
+
+    def __init__(
+        self,
+        condition: Callable[[Circuit, Circuit], bool],
+        pass_list: BasePass | Sequence[BasePass],
+    ) -> None:
+        """
+        Construct a DoThenDecide.
+
+        Args:
+            condition (Callable[[Circuit, Circuit], bool]): The condition
+                function that determines if the new circuit
+                (second parameter) after running `pass_list` should replace
+                the original circuit (first parameter). If the condition
+                returns True, then replaced the original circuit with the
+                new one.
+
+            pass_list (BasePass | Sequence[BasePass]): The pass or passes
+                to execute.
+
+        Raises:
+            ValueError: If a Sequence[BasePass] is given, but it is empty.
+        """
+        if not is_sequence(pass_list) and not isinstance(pass_list, BasePass):
+            raise TypeError(
+                'Expected Pass or sequence of Passes, got %s.'
+                % type(pass_list),
+            )
+
+        if is_sequence(pass_list):
+            truth_list = [isinstance(elem, BasePass) for elem in pass_list]
+            if not all(truth_list):
+                raise TypeError(
+                    'Expected Pass or sequence of Passes, got %s.'
+                    % type(pass_list[truth_list.index(False)]),
+                )
+            if len(pass_list) == 0:
+                raise ValueError('Expected at least one pass.')
+
+        self.condition = condition
+        self.pass_list = pass_list if is_sequence(pass_list) else [pass_list]
+
+    def run(self, circuit: Circuit, data: dict[str, Any] = {}) -> None:
+        """Perform the pass's operation, see :class:`BasePass` for more."""
+        # Perform Work
+        old_circuit = circuit.copy()
+        _logger.debug('Pass list executing...')
+        for bqskit_pass in self.pass_list:
+            bqskit_pass.run(circuit, data)
+
+        if self.condition(old_circuit, circuit):
+            _logger.debug('Accepted circuit.')
+
+        else:
+            circuit.become(old_circuit)
+            _logger.debug('Rejected circuit.')
diff --git a/bqskit/passes/control/foreach.py b/bqskit/passes/control/foreach.py
index c2858a315..24fe81ba3 100644
--- a/bqskit/passes/control/foreach.py
+++ b/bqskit/passes/control/foreach.py
@@ -7,7 +7,9 @@
 from typing import Callable
 from typing import Sequence
 
-from distributed import worker_client
+from distributed import get_client
+from distributed import rejoin
+from distributed import secede
 
 from bqskit.compiler.basepass import BasePass
 from bqskit.compiler.machine import MachineModel
@@ -157,23 +159,29 @@ def run(self, circuit: Circuit, data: dict[str, Any] = {}) -> None:
 
         # Perform Work
         if 'executor' in data:  # In Parallel
-            with worker_client() as client:
-                completed_subcircuits = []
-                completed_block_datas = []
-                futures = []
-                for subcircuit, block_data in zip(subcircuits, block_datas):
-                    future = client.submit(
-                        _sub_do_work,
-                        self.loop_body,
-                        subcircuit,
-                        block_data,
-                    )
-                    futures.append(future)
-
-                for future in futures:
-                    x, y = future.result()
-                    completed_subcircuits.append(x)
-                    completed_block_datas.append(y)
+            for block_data in block_datas:
+                block_data['executor'] = data['executor']
+            client = get_client()
+            completed_subcircuits = []
+            completed_block_datas = []
+            futures = []
+            subc_futures = client.scatter(subcircuits)
+            data_futures = client.scatter(block_datas)
+            for subcircuit, block_data in zip(subc_futures, data_futures):
+                future = client.submit(
+                    _sub_do_work,
+                    self.loop_body,
+                    subcircuit,
+                    block_data,
+                )
+                futures.append(future)
+            secede()
+            client.gather(futures)
+            rejoin()
+            for future in futures:
+                x, y = future.result()
+                completed_subcircuits.append(x)
+                completed_block_datas.append(y)
 
         else:  # Sequentially
             completed_subcircuits = []
@@ -190,6 +198,7 @@ def run(self, circuit: Circuit, data: dict[str, Any] = {}) -> None:
             # Record Data Part 2
             block_data['subcircuit_post'] = subcircuit.copy()
             block_data['loop_body_data'] = subdata
+            block_data['point'] = CircuitPoint(cycle, op.location[0])
 
             # Calculate Errors
             if self.calculate_error_bound:
diff --git a/bqskit/passes/control/paralleldo.py b/bqskit/passes/control/paralleldo.py
new file mode 100644
index 000000000..6ccbf08f0
--- /dev/null
+++ b/bqskit/passes/control/paralleldo.py
@@ -0,0 +1,125 @@
+"""This module implements the ParallelDo class."""
+from __future__ import annotations
+
+import logging
+from typing import Any
+from typing import Callable
+from typing import Sequence
+
+from distributed import as_completed
+from distributed import worker_client
+
+from bqskit.compiler.basepass import BasePass
+from bqskit.ir.circuit import Circuit
+from bqskit.utils.typing import is_sequence
+
+_logger = logging.getLogger(__name__)
+
+
+class ParallelDo(BasePass):
+    """
+    The ParallelDo class.
+
+    This is a control pass that executes a group of pass sequences in parallel.
+    """
+
+    def __init__(
+        self,
+        pass_sequences: Sequence[BasePass | Sequence[BasePass]],
+        less_than: Callable[[Circuit, Circuit], bool],
+        pick_fisrt: bool = False,
+    ) -> None:
+        """
+        Construct a ParallelDo.
+
+        Args:
+            pass_sequences (Sequence[BasePass | Sequence[BasePass]]):
+                The group of pass sequences to run in parallel.
+
+            less_than (Callable[[Circuit, Circuit], bool]): Return True
+                if the first circuit is preferred to the second one.
+                This will be used to determine which output circuit to
+                select.
+
+            pick_first (bool): If true, then the pass will complete as
+                soon as one of the subsequences finishes and will return
+                the first result.
+
+        Raises:
+            ValueError: If a Sequence[BasePass] is given, but it is empty.
+        """
+        if not is_sequence(pass_sequences):
+            raise TypeError(
+                'Expected sequence of sequences of passes, got %s.'
+                % type(pass_sequences),
+            )
+
+        self.pass_seqs = []
+        for ps in pass_sequences:
+            if not is_sequence(ps) and not isinstance(ps, BasePass):
+                raise TypeError(
+                    'Expected Pass or sequence of Passes, got %s.'
+                    % type(ps),
+                )
+
+            if is_sequence(ps):
+                truth_list = [isinstance(elem, BasePass) for elem in ps]
+                if not all(truth_list):
+                    raise TypeError(
+                        'Expected Pass or sequence of Passes, got %s.'
+                        % type(ps[truth_list.index(False)]),
+                    )
+                if len(ps) == 0:
+                    raise ValueError('Expected at least one pass.')
+
+            self.pass_seqs.append(ps if is_sequence(ps) else [ps])
+
+        self.less_than = less_than
+        self.pick_first = pick_fisrt
+
+    def run(self, circuit: Circuit, data: dict[str, Any] = {}) -> None:
+        """Perform the pass's operation, see :class:`BasePass` for more."""
+        _logger.debug('Running pass sequences in parallel.')
+
+        if 'executor' not in data:
+            raise RuntimeError(
+                'Parallel pass cannot be executed without a compiler.',
+            )
+
+        with worker_client() as client:
+            futures = []
+            for pass_seq in self.pass_seqs:
+                future = client.submit(_sub_do_work, pass_seq, circuit, data)
+                futures.append(future)
+
+            if self.pick_first:
+                _, result = next(as_completed(futures, with_results=True))
+                for future in futures:
+                    if not future.done:
+                        client.cancel(future)
+
+                circuit.become(result[0])
+                data.update(result[1])
+
+            else:
+                best_circ = None
+                best_data = None
+                for future in futures:
+                    _circ, _data = future.result()
+                    if best_circ is None or self.less_than(_circ, best_circ):
+                        best_circ = _circ
+                        best_data = _data
+                if best_circ is None or best_data is None:
+                    raise RuntimeError('No valid circuit found.')
+                circuit.become(best_circ)
+                data.update(best_data)
+
+
+def _sub_do_work(
+    loop_body: Sequence[BasePass],
+    circuit: Circuit,
+    data: dict[str, Any],
+) -> tuple[Circuit, dict[str, Any]]:
+    for loop_pass in loop_body:
+        loop_pass.run(circuit, data)
+    return circuit, data
diff --git a/bqskit/passes/partitioning/scan.py b/bqskit/passes/partitioning/scan.py
index bec24a4c5..c7486502e 100644
--- a/bqskit/passes/partitioning/scan.py
+++ b/bqskit/passes/partitioning/scan.py
@@ -3,17 +3,30 @@
 
 import logging
 from typing import Any
+from typing import Callable
+from typing import Iterator
+from typing import Sequence
+from typing import Tuple
 
 from bqskit.compiler.basepass import BasePass
 from bqskit.compiler.machine import MachineModel
 from bqskit.ir.circuit import Circuit
 from bqskit.ir.gates.circuitgate import CircuitGate
+from bqskit.ir.operation import Operation
 from bqskit.ir.region import CircuitRegion
 from bqskit.utils.typing import is_integer
 
 _logger = logging.getLogger(__name__)
 
 
+def default_scoring_fn(ops: list[Operation]) -> float:
+    """Default parition scoring function based on QGo."""
+    score = 0.0
+    for op in ops:
+        score += (op.num_qudits - 1) * 100 + 1
+    return score
+
+
 class ScanPartitioner(BasePass):
     """
     The ScanPartitioner Pass.
@@ -28,9 +41,8 @@ class ScanPartitioner(BasePass):
 
     def __init__(
         self,
-        block_size: int = 3,
-        single_gate_score: int = 1,
-        multi_gate_score: int = 1000,  # TODO: Pass callable scoring_fn instead
+        block_size: int,
+        scoring_fn: Callable[[list[Operation]], float] = default_scoring_fn,
     ) -> None:
         """
         Construct a ScanPartitioner.
@@ -39,13 +51,10 @@ def __init__(
             block_size (int): Maximum size of partitioned blocks.
                 (Default: 3)
 
-            single_gate_score (int): When evaluating potential blocks,
-                use this number to score included single-qudit gates.
-                (Default: 1)
-
-            multi_gate_score (int): When evaluating potential blocks,
-                use this number to score included multi-qudit gates.
-                (Default: 1000)
+            scoring_fn (Callable[[list[Operation]], float]): This function
+                is used to score potential blocks. This takes in a list of
+                operations in the potential block and returns a float.
+                (Default: default_scoring_fn)
 
         Raises:
             ValueError: If `block_size` is less than 2.
@@ -61,21 +70,8 @@ def __init__(
                 f'Expected block_size to be greater than 2, got {block_size}.',
             )
 
-        if not is_integer(single_gate_score):
-            raise TypeError(
-                'Expected integer for single_gate_score, '
-                f'got {type(single_gate_score)}.',
-            )
-
-        if not is_integer(multi_gate_score):
-            raise TypeError(
-                'Expected integer for multi_gate_score, '
-                f'got {type(multi_gate_score)}.',
-            )
-
         self.block_size = block_size
-        self.single_gate_score = single_gate_score
-        self.multi_gate_score = multi_gate_score
+        self.scoring_fn = scoring_fn
 
     def run(self, circuit: Circuit, data: dict[str, Any] = {}) -> None:
         """Perform the pass's operation, see :class:`BasePass` for more."""
@@ -91,166 +87,297 @@ def run(self, circuit: Circuit, data: dict[str, Any] = {}) -> None:
             })
             return
 
-        # If a MachineModel is provided in the data dict, it will be used.
-        # Otherwise all-to-all connectivity is assumed.
-        model = None
+        # Cache maximum cycles in circuit
+        num_cycles = circuit.num_cycles
 
-        if 'machine_model' in data:
-            model = data['machine_model']
+        # All groups of qudits that can exist in viable partitions
+        qudit_groups = self.calculate_qudit_groups(circuit)
 
-        if (
-            not isinstance(model, MachineModel)
-            or model.num_qudits < circuit.num_qudits
-        ):
-            _logger.warning(
-                'MachineModel not specified or invalid;'
-                ' defaulting to all-to-all.',
-            )
-            model = MachineModel(circuit.num_qudits)
-
-        # Find all connected, `block_size`-sized groups of qudits
-        # NOTE: This assumes circuit and topology qudit numbers are equal
-        qudit_groups = model.get_locations(self.block_size)
-        # Prune unused qudit groups
-        used_qudits = [
-            q for q in range(circuit.num_qudits)
-            if not circuit.is_qudit_idle(q)
-        ]
-        for qudit_group in qudit_groups:
-            if all([q not in used_qudits for q in qudit_group]):
-                qudit_groups.remove(qudit_group)
+        # Map from qudit indices to a list of groups they exist in
+        qudit_group_map = self.calculate_qudit_group_map(qudit_groups)
+
+        # Map from qudit groups to current best block for that group
+        potential_blocks = self.calculate_initial_blocks(qudit_groups, circuit)
 
         # divider splits the circuit into partitioned and unpartitioned spaces.
-        active_qudits = circuit.active_qudits
-        num_cycles = circuit.num_cycles
         divider = [
-            0 if q in active_qudits else num_cycles
+            0 if q in qudit_group_map.keys() else num_cycles
             for q in range(circuit.num_qudits)
         ]
 
-        # Form regions until there are no more gates to partition
+        # Stores the selected blocks as regions
         regions: list[CircuitRegion] = []
+
+        # Form regions until there are no more gates to partition
         while any(cycle < num_cycles for cycle in divider):
 
-            # Move past/skip any gates that are larger than block size
-            qudits_to_increment: list[int] = []
-            for qudit, cycle in enumerate(divider):
-                if qudit in qudits_to_increment or cycle >= num_cycles:
-                    continue
-
-                if not circuit.is_point_idle((cycle, qudit)):
-                    op = circuit[cycle, qudit]
-                    if len(op.location) > self.block_size:
-                        if all(divider[q] == cycle for q in op.location):
-                            qudits_to_increment.extend(op.location)
-                            _logger.warning(
-                                'Skipping gate larger than block size.',
-                            )
-
-            if len(qudits_to_increment) > 0:
-                regions.append(
-                    CircuitRegion(
-                        {
-                            qudit: (divider[qudit], divider[qudit])
-                            for qudit in qudits_to_increment
-                        },
-                    ),
-                )
+            # Select best block from current potential blocks
+            best_region = self.find_best_block(potential_blocks)
+            regions.append(best_region)
+            _logger.info(f'Just formed region: {best_region}')
 
-            for qudit in qudits_to_increment:
-                divider[qudit] += 1
-
-            # Skip any idle qudit-cycles
-            amount_to_add_to_each_qudit = [0 for _ in range(circuit.num_qudits)]
-            for qudit, cycle in enumerate(divider):
-                while (
-                    cycle < num_cycles
-                    and circuit.is_point_idle((cycle, qudit))
-                ):
-                    amount_to_add_to_each_qudit[qudit] += 1
-                    cycle += 1
-
-            for qudit, amount in enumerate(amount_to_add_to_each_qudit):
-                divider[qudit] += amount
-
-            # Find the scores of the qudit groups.
-            best_score = None
-            best_region = None
-
-            for qudit_group in qudit_groups:
-
-                ops_and_cycles = circuit.operations_with_cycles(
-                    qudits_or_region=CircuitRegion({
-                        qudit_index: (divider[qudit_index], num_cycles)
-                        for qudit_index in qudit_group
-                    }),
+            # Update divider
+            for qudit, interval in best_region.items():
+                divider[qudit] = interval[1] + 1
+
+            # Update adjacent groups
+            adjacent_groups = set()
+            for qudit in best_region.keys():
+                adjacent_groups.update(qudit_group_map[qudit])
+
+            for qudit_group in adjacent_groups:
+                starting_cycles = [divider[q] for q in qudit_group]
+                new_block = self.calculate_block(
+                    qudit_group, circuit, starting_cycles,
                 )
+                potential_blocks[qudit_group] = new_block
+
+        # Form the partitioned circuit
+        circuit.become(self.fold_circuit(circuit, regions))
 
-                in_qudits = list(q for q in qudit_group)
-                stopped_cycles = {q: num_cycles for q in qudit_group}
-                score = 0
-
-                for cycle, op in ops_and_cycles:
-                    if len(op.location.union(in_qudits)) != len(in_qudits):
-                        for qudit_index in op.location.intersection(in_qudits):
-                            stopped_cycles[qudit_index] = cycle
-                            in_qudits.remove(qudit_index)
-                    else:
-                        if len(op.location) > 1:
-                            score += self.multi_gate_score
-                        else:
-                            score += self.single_gate_score
-
-                    if len(in_qudits) == 0:
+    def fold_circuit(
+        self,
+        circuit: Circuit,
+        regions: list[CircuitRegion],
+    ) -> Circuit:
+        """Partition `circuit` into blocks described in `regions`."""
+        folded_circuit = Circuit(circuit.num_qudits, circuit.radixes)
+
+        for region in regions:
+            ops_and_cycles = list({
+                (point[0], circuit[point])
+                for point in region.points
+                if not circuit.is_point_idle(point)
+            })
+            ops_and_cycles.sort(key=lambda x: (x[0], *x[1].location))
+            ops = [op for _, op in ops_and_cycles]
+            qudits = list(set(sum((tuple(op.location) for op in ops), ())))
+            qudits = sorted(qudits)
+            radixes = [circuit.radixes[q] for q in qudits]
+
+            cgc = Circuit(len(radixes), radixes)
+            for op in ops:
+                location = [qudits.index(q) for q in op.location]
+                cgc.append(Operation(op.gate, location, op.params))
+
+            folded_circuit.append_gate(
+                CircuitGate(cgc, True),
+                qudits,
+                cgc.params,
+            )
+
+        return folded_circuit
+
+    def calculate_qudit_groups(self, circuit: Circuit) -> list[tuple[int, ...]]:
+        """Calculate the groups of qudits that grouped in a partition."""
+        # Get initial set from circuit connectivity
+        model = MachineModel(circuit.num_qudits, circuit.coupling_graph)
+        qudit_groups: list[tuple[int, ...]] = []
+        for i in reversed(range(self.block_size)):
+            potential_groups = [tuple(x) for x in model.get_locations(i + 1)]
+            for potential_group in potential_groups:
+                should_add_group = True
+                for qudit_group in qudit_groups:
+                    if all(x in qudit_group for x in potential_group):
+                        should_add_group = False
                         break
+                if should_add_group:
+                    qudit_groups.append(potential_group)
 
-                if best_score is None or score > best_score:
-                    best_score = score
-                    best_region = CircuitRegion({
-                        qudit: (
-                            divider[qudit],
-                            stopped_cycles[qudit] - 1,
-                        )
-                        for qudit in qudit_group
-                        if stopped_cycles[qudit] - 1 >= divider[qudit]
-                    })
-
-            if best_score is None or best_region is None:
-                raise RuntimeError('No valid block found.')
-
-            _logger.info('Found block with score: %d.' % (best_score))
-            regions.append(best_region)
+        # Prune groups
+        active_qudits = circuit.active_qudits
+        qudits_in_groups: set[int] = set()
+        qudit_groups_to_remove = []
+        qudit_groups_to_append = []
+        for qudit_group in qudit_groups:
+            if all(q in active_qudits for q in qudit_group):
+                qudits_in_groups.update(qudit_group)
+            else:
+                qudit_groups_to_remove.append(qudit_group)
+                group = tuple(q for q in qudit_group if q in active_qudits)
+                if len(group) > 0:
+                    qudit_groups_to_append.append(group)
+                    qudits_in_groups.update(group)
 
-            # Update divider
-            for qudit_index in best_region:
-                divider[qudit_index] = best_region[qudit_index].upper + 1
+        for qudit_group in qudit_groups_to_remove:
+            qudit_groups.remove(qudit_group)
 
-        # Fold the circuit
-        folded_circuit = Circuit(circuit.num_qudits, circuit.radixes)
-        # Option to keep a block's idle qudits as part of the CircuitGate
-        if 'keep_idle_qudits' in data and data['keep_idle_qudits'] is True:
-            for region in regions:
-                small_region = circuit.downsize_region(region)
-                cgc = circuit.get_slice(small_region.points)
-                if len(region.location) > len(small_region.location):
-                    for i in range(len(region.location)):
-                        if region.location[i] not in small_region.location:
-                            cgc.insert_qudit(i)
-                folded_circuit.append_gate(
-                    CircuitGate(cgc, True),
-                    sorted(list(region.keys())),
-                    list(cgc.params),
-                )
-        else:
-            for region in regions:
-                region = circuit.downsize_region(region)
-                if 0 < len(region) <= self.block_size:
-                    cgc = circuit.get_slice(region.points)
-                    folded_circuit.append_gate(
-                        CircuitGate(cgc, True),
-                        sorted(list(region.keys())),
-                        list(cgc.params),
+        for qudit_group in qudit_groups_to_append:
+            qudit_groups.append(qudit_group)
+
+        # Add groups for individual active qudits
+        for q in active_qudits:
+            if q not in qudits_in_groups:
+                qudit_groups.append((q,))
+
+        return qudit_groups
+
+    def calculate_qudit_group_map(
+        self,
+        qudit_groups: list[tuple[int, ...]],
+    ) -> dict[int, list[tuple[int, ...]]]:
+        """Calculate a map from qudits to qudit_groups containing that qudit."""
+        qudit_to_group_map: dict[int, list[tuple[int, ...]]] = {}
+        for qudit_group in qudit_groups:
+            for qudit in qudit_group:
+                if qudit not in qudit_to_group_map:
+                    qudit_to_group_map[qudit] = []
+                qudit_to_group_map[qudit].append(qudit_group)
+        return qudit_to_group_map
+
+    def calculate_initial_blocks(
+        self,
+        qudit_groups: Sequence[tuple[int, ...]],
+        circuit: Circuit,
+    ) -> dict[tuple[int, ...], tuple[CircuitRegion, list[Operation]]]:
+        """Calculate the initial blocks for all qudit groups."""
+
+        blocks: dict[
+            tuple[int, ...],
+            tuple[CircuitRegion, list[Operation]],
+        ] = {}
+
+        for qg in qudit_groups:
+            blocks[qg] = self.calculate_block(qg, circuit, [0] * len(qg))
+
+        return blocks
+
+    def find_best_block(
+        self,
+        potential_blocks: dict[
+            tuple[int, ...],
+            tuple[CircuitRegion, list[Operation]],
+        ],
+    ) -> CircuitRegion:
+        """Find and return the current best scoring block."""
+        return sorted(
+            list(potential_blocks.values()),
+            key=lambda x: self.scoring_fn(x[1]),
+        )[-1][0]
+
+    def calculate_block(
+        self,
+        qudit_group: Sequence[int],
+        circuit: Circuit,
+        starting_cycles: Sequence[int],
+    ) -> tuple[CircuitRegion, list[Operation]]:
+        """
+        Calculate the best block for `qudit_group` right of the divider.
+
+        Args:
+            qudit_group (Sequence[int]): The block's qudits.
+
+            circuit (Circuit): The circuit to form a block in.
+
+            starting_cycles (Sequence[int]): Where to start scanning
+                right from.
+
+        Returns:
+            (tuple[CircuitRegion, list[Operation]]): The formed region
+                and the operations in the block.
+        """
+        in_qudits = list(q for q in qudit_group)
+        stopped_cycles = {q: circuit.num_cycles for q in qudit_group}
+        op_list: list[Operation] = []
+
+        iter = self.FastRegionIterator(qudit_group, starting_cycles, circuit)
+        for cycle, op in iter:
+            if any(op_q not in in_qudits for op_q in op.location):
+                if op.num_qudits > self.block_size:
+                    raise RuntimeError(
+                        'ScanPartitioner cannot handle gates larger than'
+                        ' block size. You may want to use the '
+                        'QuickPartitioner.',
                     )
-                else:
-                    folded_circuit.extend(circuit[region])
-        circuit.become(folded_circuit)
+
+                for op_q in op.location:
+                    if op_q in in_qudits:
+                        stopped_cycles[op_q] = cycle
+                        in_qudits.remove(op_q)
+            else:
+                op_list.append(op)
+
+            if len(in_qudits) == 0:
+                break
+
+        return (
+            CircuitRegion({
+                q: (starting_cycles[i], stopped_cycles[q] - 1)
+                for i, q in enumerate(qudit_group)
+                if stopped_cycles[q] - 1 >= starting_cycles[i]
+            }),
+            op_list,
+        )
+
+    class FastRegionIterator(Iterator[Tuple[int, Operation]]):
+        """A circuit iterator designed to be efficient for the ScanPartitioner's
+        use case."""
+
+        def __init__(
+            self,
+            qudits: Sequence[int],
+            starting_cycles: Sequence[int],
+            circuit: Circuit,
+        ) -> None:
+            """Construct a FastRegionIterator."""
+            self.qudits = qudits
+            self.starting_cycles = {
+                q: s for q, s in zip(
+                    qudits, starting_cycles,
+                )
+            }
+            self.circuit = circuit
+            self.num_cycles = circuit.num_cycles
+
+            self.inactive = list(
+                sorted(zip(qudits, starting_cycles), key=lambda x: x[1]),
+            )
+            self.active = []
+            min_cycle = self.inactive[0][1]
+            for qudit, starting_cycle in self.inactive:
+                if starting_cycle == min_cycle:
+                    self.active.append(qudit)
+            for qudit in self.active:
+                self.inactive.pop(0)
+
+            self.cycle = self.starting_cycles[self.active[0]]
+            self.qudit_index = 0
+            self.qudit = self.active[self.qudit_index]
+            self.qudits_to_skip: set[int] = set()
+
+        def __iter__(self) -> Iterator[tuple[int, Operation]]:
+            return self
+
+        def step(self) -> None:
+            """Move the iterator one forward."""
+            self.qudit_index += 1
+
+            if self.qudit_index >= len(self.active):
+                self.qudit_index = 0
+                self.cycle += 1
+                self.qudits_to_skip.clear()
+                while True:
+                    if len(self.inactive) == 0:
+                        break
+                    if self.inactive[0][1] != self.cycle:
+                        break
+                    self.active.append(self.inactive[0][0])
+                    self.inactive.pop(0)
+
+            self.qudit = self.active[self.qudit_index]
+
+        def __next__(self) -> tuple[int, Operation]:
+            """Get the next operation."""
+            if self.cycle >= self.num_cycles:
+                raise StopIteration
+            op = self.circuit._circuit[self.cycle][self.qudit]
+            while op is None:
+                self.step()
+                if self.cycle >= self.num_cycles:
+                    raise StopIteration
+                op = self.circuit._circuit[self.cycle][self.qudit]
+
+            self.qudits_to_skip.update(op.location)
+            cycle_to_return = self.cycle
+            self.step()
+            return cycle_to_return, op
diff --git a/bqskit/passes/processing/__init__.py b/bqskit/passes/processing/__init__.py
index 2b09ca101..1d1e8170d 100644
--- a/bqskit/passes/processing/__init__.py
+++ b/bqskit/passes/processing/__init__.py
@@ -2,6 +2,11 @@
 from __future__ import annotations
 
 from bqskit.passes.processing.scan import ScanningGateRemovalPass
+from bqskit.passes.processing.substitute import SubstitutePass
 from bqskit.passes.processing.window import WindowOptimizationPass
 
-__all__ = ['ScanningGateRemovalPass', 'WindowOptimizationPass']
+__all__ = [
+    'ScanningGateRemovalPass',
+    'SubstitutePass',
+    'WindowOptimizationPass',
+]
diff --git a/bqskit/passes/processing/substitute.py b/bqskit/passes/processing/substitute.py
new file mode 100644
index 000000000..747e978f3
--- /dev/null
+++ b/bqskit/passes/processing/substitute.py
@@ -0,0 +1,139 @@
+"""This module implements the SubstitutePass."""
+from __future__ import annotations
+
+import itertools as it
+import logging
+from typing import Any
+from typing import Callable
+
+from bqskit.compiler.basepass import BasePass
+from bqskit.ir.circuit import Circuit
+from bqskit.ir.gate import Gate
+from bqskit.ir.operation import Operation
+from bqskit.ir.opt.cost.functions import HilbertSchmidtResidualsGenerator
+from bqskit.ir.opt.cost.generator import CostFunctionGenerator
+from bqskit.utils.typing import is_real_number
+_logger = logging.getLogger(__name__)
+
+
+class SubstitutePass(BasePass):
+    """
+    The SubstitutePass class.
+
+    The substitute pass will attempt to use instantiation to replace gates with
+    other gates.
+    """
+
+    def __init__(
+        self,
+        collection_filter: Callable[[Operation], bool],
+        gate: Gate,
+        success_threshold: float = 1e-10,
+        cost: CostFunctionGenerator = HilbertSchmidtResidualsGenerator(),
+        instantiate_options: dict[str, Any] = {},
+    ) -> None:
+        """
+        Construct a SubstitutePass.
+
+        Args:
+            collection_filter (Callable[[Operation], bool]):
+                A predicate that determines which operations should have
+                substitution attempted on them.
+
+            gate (Gate): The gate to try substitute in.
+
+            success_threshold (float): The distance threshold that
+                determines successful termintation. Measured in cost
+                described by the hilbert schmidt cost function.
+                (Default: 1e-10)
+
+            cost (CostFunction | None): The cost function that determines
+                successful removal of a gate.
+                (Default: HilbertSchmidtResidualsGenerator())
+
+            instantiate_options (dict[str: Any]): Options passed directly
+                to circuit.instantiate when instantiating circuit
+                templates. (Default: {})
+        """
+        if not callable(collection_filter):
+            raise TypeError(
+                'Expected callable method that maps Operations to booleans for'
+                ' collection_filter, got %s.' % type(self.collection_filter),
+            )
+
+        if not isinstance(gate, Gate):
+            raise TypeError(f'Expected a gate, got {type(gate)}')
+
+        if not is_real_number(success_threshold):
+            raise TypeError(
+                'Expected real number for success_threshold'
+                ', got %s' % type(success_threshold),
+            )
+
+        if not isinstance(cost, CostFunctionGenerator):
+            raise TypeError(
+                'Expected cost to be a CostFunctionGenerator, got %s'
+                % type(cost),
+            )
+
+        if not isinstance(instantiate_options, dict):
+            raise TypeError(
+                'Expected dictionary for instantiate_options, got %s.'
+                % type(instantiate_options),
+            )
+
+        self.collection_filter = collection_filter
+        self.gate = gate
+        self.success_threshold = success_threshold
+        self.cost = cost
+        self.instantiate_options: dict[str, Any] = {
+            'dist_tol': self.success_threshold,
+        }
+        self.instantiate_options.update(instantiate_options)
+
+    def run(self, circuit: Circuit, data: dict[str, Any] = {}) -> None:
+        """Perform the pass's operation, see :class:`BasePass` for more."""
+
+        # Update the instantiater to run in parallel
+        instantiate_options = self.instantiate_options
+        if 'executor' in data:
+            instantiate_options['parallel'] = True
+
+        # Collect locations in circuit where self.init_gate exists
+        points = [
+            (cycle, op.location[0])
+            for cycle, op in circuit.operations_with_cycles()
+            if self.collection_filter(op)
+        ]
+
+        # Check substitutions are valid
+        for point in points:
+            if circuit[point].num_qudits < self.gate.num_qudits:
+                raise RuntimeError(
+                    f'Cannot substitute {circuit[point]}'
+                    f' with {self.gate}.',
+                )
+
+        # Attempt substitution
+        target = circuit.get_unitary()
+        shift = 0
+        for point in points:
+            point = (point[0] - shift, point[1])
+            _logger.info(f'Attempting substitution for operation at {point}.')
+            _logger.info(f'Operation: {circuit[point]}')
+
+            qudits = circuit[point].location
+            locs = it.combinations(qudits, self.gate.num_qudits)
+            og_cycle_count = circuit.num_cycles
+
+            for loc in locs:
+                _logger.debug(f'Trying location: {loc}')
+                circuit_copy = circuit.copy()
+                circuit_copy.replace_gate(point, self.gate, loc)
+                circuit_copy.instantiate(target, **self.instantiate_options)
+
+                if self.cost(circuit_copy, target) < self.success_threshold:
+                    _logger.info('Successfully substituted operation.')
+                    circuit.become(circuit_copy)
+                    shift += og_cycle_count - circuit.num_cycles
+                    break
diff --git a/bqskit/passes/search/generators/__init__.py b/bqskit/passes/search/generators/__init__.py
index f3410c1bf..d96eda3f6 100644
--- a/bqskit/passes/search/generators/__init__.py
+++ b/bqskit/passes/search/generators/__init__.py
@@ -3,5 +3,11 @@
 
 from bqskit.passes.search.generators.seed import SeedLayerGenerator
 from bqskit.passes.search.generators.simple import SimpleLayerGenerator
+from bqskit.passes.search.generators.stair import StairLayerGenerator
 
-__all__ = ['SimpleLayerGenerator', 'SeedLayerGenerator']
+__all__ = [
+    'SimpleLayerGenerator',
+    'SeedLayerGenerator',
+    'StairLayerGenerator',
+    'MiddleOutLayerGenerator',
+]
diff --git a/bqskit/passes/search/generators/middleout.py b/bqskit/passes/search/generators/middleout.py
new file mode 100644
index 000000000..e65466093
--- /dev/null
+++ b/bqskit/passes/search/generators/middleout.py
@@ -0,0 +1,75 @@
+"""This module implements the MiddleOutLayerGenerator class."""
+from __future__ import annotations
+
+import itertools as it
+import logging
+from typing import Any
+
+from bqskit.compiler.basepass import BasePass
+from bqskit.ir.circuit import Circuit
+from bqskit.passes.search.generators.seed import SeedLayerGenerator
+from bqskit.passes.search.generators.simple import SimpleLayerGenerator
+
+_logger = logging.getLogger(__name__)
+
+
+# TODO: choose a better name?
+class MiddleOutLayerGenerator(SimpleLayerGenerator):
+    """Layer Generator for search that adds gates at each time step."""
+    generated_circuits: set[int]
+
+    def __init__(self) -> None:
+        """Construct a MiddleOutLayerGenerator."""
+        self.generated_circuits = set()
+        super().__init__()
+
+    def gen_successors(
+        self,
+        circuit: Circuit,
+        data: dict[str, Any],
+    ) -> list[Circuit]:
+        """
+        Generate the successors of a circuit node.
+
+        Raises:
+            ValueError: If circuit is a single-qudit circuit.
+        """
+
+        if not isinstance(circuit, Circuit):
+            raise TypeError('Expected circuit, got %s.' % type(circuit))
+
+        if circuit.num_qudits < 2:
+            raise ValueError('Cannot expand a single-qudit circuit.')
+
+        # Get the machine model
+        model = BasePass.get_model(circuit, data)
+
+        # Generate successors
+        successors = []
+        for (edge, cycle) in it.product(
+                model.coupling_graph, range(circuit.num_cycles),
+        ):
+            successor = circuit.copy()
+            successor.insert_gate(
+                cycle, self.two_qudit_gate, [
+                    edge[0], edge[1],
+                ],
+            )
+            successor.insert_gate(cycle, self.single_qudit_gate_1, edge[0])
+            successor.insert_gate(cycle, self.single_qudit_gate_2, edge[1])
+            # TODO: evaluate whether or not this hash method is efficient
+            hash = SeedLayerGenerator.hash_structure(successor)
+            if hash not in self.generated_circuits:
+                self.generated_circuits.add(hash)
+                successors.append(successor)
+
+        for edge in model.coupling_graph:
+            successor = circuit.copy()
+            successor.append_gate(self.two_qudit_gate, [edge[0], edge[1]])
+            successor.append_gate(self.single_qudit_gate_1, edge[0])
+            successor.append_gate(self.single_qudit_gate_2, edge[1])
+            if hash not in self.generated_circuits:
+                self.generated_circuits.add(hash)
+                successors.append(successor)
+
+        return successors
diff --git a/bqskit/passes/search/generators/seed.py b/bqskit/passes/search/generators/seed.py
index b2b5a5526..bbc4b99cf 100644
--- a/bqskit/passes/search/generators/seed.py
+++ b/bqskit/passes/search/generators/seed.py
@@ -23,7 +23,7 @@ def __init__(
         num_removed: int = 1,
     ) -> None:
         """
-        Construct a SimpleLayerGenerator.
+        Construct a SeedLayerGenerator.
 
         Args:
             seed (Circuit): The seed to start from.
diff --git a/bqskit/passes/search/generators/stair.py b/bqskit/passes/search/generators/stair.py
new file mode 100644
index 000000000..e7aa66fdd
--- /dev/null
+++ b/bqskit/passes/search/generators/stair.py
@@ -0,0 +1,79 @@
+"""This module implements the StairLayerGenerator class."""
+from __future__ import annotations
+
+import logging
+from typing import Any
+
+from bqskit.ir.circuit import Circuit
+from bqskit.ir.gate import Gate
+from bqskit.passes.search.generator import LayerGenerator
+from bqskit.qis.state.state import StateVector
+from bqskit.qis.unitary.unitarymatrix import UnitaryMatrix
+
+
+_logger = logging.getLogger(__name__)
+
+
+class StairLayerGenerator(LayerGenerator):
+    """Layer Generator for search that builds circuits from a single gate."""
+
+    def __init__(self, gate: Gate) -> None:
+        """
+        Construct a StairLayerGenerator.
+
+        Args:
+            gate (Gate): The gate to build from.
+        """
+        if not isinstance(gate, Gate):
+            raise TypeError(f'Expected Gate for gate, got {type(gate)}')
+
+        self.gate = gate
+
+    def gen_initial_layer(
+        self,
+        target: UnitaryMatrix | StateVector,
+        data: dict[str, Any],
+    ) -> Circuit:
+        """
+        Generate the initial layer, see LayerGenerator for more.
+
+        Raises:
+            ValueError: If `target` has a size or radix mismatch with
+                `self.seed`.
+        """
+
+        if not isinstance(target, (UnitaryMatrix, StateVector)):
+            raise TypeError(
+                'Expected unitary or state, got %s.' % type(target),
+            )
+
+        init_circuit = Circuit(target.num_qudits, target.radixes)
+
+        # Place RXX Gates on consective pairs of qudits
+        nq = self.gate.num_qudits
+        for i in range(init_circuit.num_qudits - (nq - 1)):
+            init_circuit.append_gate(self.gate, [i + j for j in range(nq)])
+
+        return init_circuit
+
+    def gen_successors(
+        self,
+        circuit: Circuit,
+        data: dict[str, Any],
+    ) -> list[Circuit]:
+        """
+        Generate the successors of a circuit node.
+
+        Raises:
+            ValueError: If circuit is a single-qudit circuit.
+        """
+
+        if not isinstance(circuit, Circuit):
+            raise TypeError('Expected circuit, got %s.' % type(circuit))
+
+        nq = self.gate.num_qudits
+        successor = circuit.copy()
+        for i in range(circuit.num_qudits - (nq - 1)):
+            successor.append_gate(self.gate, [i + j for j in range(nq)])
+
+        return [successor]
diff --git a/bqskit/passes/synthesis/__init__.py b/bqskit/passes/synthesis/__init__.py
index a78c31194..0cc3ebaa3 100644
--- a/bqskit/passes/synthesis/__init__.py
+++ b/bqskit/passes/synthesis/__init__.py
@@ -2,7 +2,6 @@
 from __future__ import annotations
 
 from bqskit.passes.synthesis.leap import LEAPSynthesisPass
-from bqskit.passes.synthesis.leap import OptimizedLEAPPass
 from bqskit.passes.synthesis.qfast import QFASTDecompositionPass
 from bqskit.passes.synthesis.qpredict import QPredictDecompositionPass
 from bqskit.passes.synthesis.qsearch import QSearchSynthesisPass
@@ -10,7 +9,6 @@
 
 __all__ = [
     'LEAPSynthesisPass',
-    'OptimizedLEAPPass',
     'QFASTDecompositionPass',
     'QPredictDecompositionPass',
     'QSearchSynthesisPass',
diff --git a/bqskit/passes/synthesis/leap.py b/bqskit/passes/synthesis/leap.py
index 3cfcc0873..5952a9e54 100644
--- a/bqskit/passes/synthesis/leap.py
+++ b/bqskit/passes/synthesis/leap.py
@@ -5,17 +5,10 @@
 from typing import Any
 
 import numpy as np
-from dask.distributed import as_completed
-from dask.distributed import worker_client
-from qsearch import assemblers
-from qsearch import leap_compiler
-from qsearch import options
-from qsearch import parallelizers
-from qsearch import post_processing
+from dask.distributed import get_client
 from scipy.stats import linregress
 
 from bqskit.ir.circuit import Circuit
-from bqskit.ir.lang.qasm2 import OPENQASM2Language
 from bqskit.ir.opt.cost.functions import HilbertSchmidtResidualsGenerator
 from bqskit.ir.opt.cost.generator import CostFunctionGenerator
 from bqskit.passes.search.frontier import Frontier
@@ -48,6 +41,8 @@ def __init__(
         success_threshold: float = 1e-10,
         cost: CostFunctionGenerator = HilbertSchmidtResidualsGenerator(),
         max_layer: int | None = None,
+        store_partial_solutions: bool = False,
+        partials_per_depth: int = 25,
         min_prefix_size: int = 3,
         instantiate_options: dict[str, Any] = {},
     ) -> None:
@@ -73,7 +68,14 @@ def __init__(
 
             max_layer (int): The maximum number of layers to append without
                 success before termination. If left as None it will default
-                 to unlimited. (Default: None)
+                to unlimited. (Default: None)
+
+            store_partial_solutions (bool): Whether to store partial solutions
+                at different depths inside of the data dict. (Default: False)
+
+            partials_per_depth (int): The maximum number of partials
+                to store per search depth. No effect if
+                `store_partial_solutions` is False. (Default: 25)
 
             min_prefix_size (int): The minimum number of layers needed
                 to prefix the circuit.
@@ -143,21 +145,21 @@ def __init__(
         self.cost = cost
         self.max_layer = max_layer
         self.min_prefix_size = min_prefix_size
-        self.instantiate_options: dict[str, Any] = {'cost_fn_gen': self.cost}
+        self.instantiate_options: dict[str, Any] = {
+            'cost_fn_gen': self.cost,
+        }
         self.instantiate_options.update(instantiate_options)
+        self.store_partial_solutions = store_partial_solutions
+        self.partials_per_depth = partials_per_depth
 
     def synthesize(self, utry: UnitaryMatrix, data: dict[str, Any]) -> Circuit:
         """Synthesize `utry`, see :class:`SynthesisPass` for more."""
         frontier = Frontier(utry, self.heuristic_function)
         data['window_markers'] = []
 
-        instantiate_options = self.instantiate_options
-        if 'executor' in data:
-            instantiate_options['parallel'] = True
-
         # Seed the search with an initial layer
         initial_layer = self.layer_gen.gen_initial_layer(utry, data)
-        initial_layer.instantiate(utry, **instantiate_options)
+        initial_layer.instantiate(utry, **self.instantiate_options)
         frontier.add(initial_layer, 0)
 
         # Track best circuit, initially the initial layer
@@ -168,6 +170,8 @@ def synthesize(self, utry: UnitaryMatrix, data: dict[str, Any]) -> Circuit:
         leap_data['best_dists'] = [leap_data['best_dist']]
         leap_data['best_layers'] = [0]
         leap_data['last_prefix_layer'] = 0
+        if self.store_partial_solutions:
+            leap_data['psols'] = {}
         _logger.info(
             'Search started, initial layer has cost: %e.' %
             leap_data['best_dist'],
@@ -177,38 +181,51 @@ def synthesize(self, utry: UnitaryMatrix, data: dict[str, Any]) -> Circuit:
             _logger.info('Successful synthesis.')
             return initial_layer
 
-        while not frontier.empty():
-            top_circuit, layer = frontier.pop()
-
-            # Generate successors and evaluate each
-            successors = self.layer_gen.gen_successors(top_circuit, data)
-
-            if 'executor' in data:  # In Parallel
-                with worker_client() as client:
-                    futures = client.map(
-                        Circuit.instantiate,
-                        successors,
-                        pure=False,
-                        target=utry,
-                        **instantiate_options,
-                    )
-                    for _, circuit in as_completed(futures, with_results=True):
-                        if self.evaluate_node(
-                            circuit,
-                            utry,
-                            data,
-                            frontier,
-                            layer,
-                            leap_data,
-                        ):
-                            for future in futures:
-                                if not future.done:
-                                    client.cancel(future)
-                            return circuit
-
-            else:  # Sequentially
+        if 'executor' in data:  # In Parallel
+            client = get_client()
+            while not frontier.empty():
+                top_circuit, layer = frontier.pop()
+
+                # Generate successors and evaluate each
+                successors = self.layer_gen.gen_successors(top_circuit, data)
+
+                # Submit instantiate jobs
+                futures = self.batched_instantiate(
+                    successors,
+                    utry,
+                    client,
+                    **self.instantiate_options,
+                )
+
+                # Wait for and gather results
+                circuits = self.gather_best_results(
+                    futures,
+                    client,
+                    self.cost.calc_cost,
+                    utry,
+                )
+
+                for circuit in circuits:
+                    if self.evaluate_node(
+                        circuit,
+                        utry,
+                        data,
+                        frontier,
+                        layer,
+                        leap_data,
+                    ):
+                        if self.store_partial_solutions:
+                            data['psols'] = leap_data['psols']
+                        return circuit
+
+        else:  # Sequentially
+            while not frontier.empty():
+                top_circuit, layer = frontier.pop()
+                # Generate successors and evaluate each
+                successors = self.layer_gen.gen_successors(top_circuit, data)
+
                 for circuit in successors:
-                    circuit.instantiate(utry, **instantiate_options)
+                    circuit.instantiate(utry, **self.instantiate_options)
                     if self.evaluate_node(
                         circuit,
                         utry,
@@ -217,6 +234,8 @@ def synthesize(self, utry: UnitaryMatrix, data: dict[str, Any]) -> Circuit:
                         layer,
                         leap_data,
                     ):
+                        if self.store_partial_solutions:
+                            data['psols'] = leap_data['psols']
                         return circuit
 
         _logger.info('Frontier emptied.')
@@ -227,6 +246,8 @@ def synthesize(self, utry: UnitaryMatrix, data: dict[str, Any]) -> Circuit:
                 else 's', leap_data['best_dist'],
             ),
         )
+        if self.store_partial_solutions:
+            data['psols'] = leap_data['psols']
 
         return leap_data['best_circ']
 
@@ -267,6 +288,16 @@ def evaluate_node(
                 if self.max_layer is None or layer + 1 < self.max_layer:
                     frontier.add(circuit, layer + 1)
 
+        if self.store_partial_solutions:
+            if layer not in leap_data['psols']:
+                leap_data['psols'][layer] = []
+
+            leap_data['psols'][layer].append((circuit.copy(), dist))
+
+            if len(leap_data['psols'][layer]) > self.partials_per_depth:
+                leap_data['psols'][layer].sort(key=lambda x: x[1])
+                del leap_data['psols'][layer][-1]
+
         if self.max_layer is None or layer + 1 < self.max_layer:
             frontier.add(circuit, layer + 1)
         return False
@@ -348,38 +379,3 @@ def check_leap_condition(
 
         layers_added = new_layer - leap_data['last_prefix_layer']
         return delta < 0 and layers_added >= self.min_prefix_size
-
-
-# TODO: Improve quality of LEAPSynthesisPass and deprecate OptimizedLEAPPass
-
-class OptimizedLEAPPass(SynthesisPass):
-    """
-    An optimized version of LEAP.
-
-    This version of the LEAP algorithm is a translation layer to the
-    original implementation of LEAP. This version has more features
-    and generally produces higher-quality results than the LEAPSynthesisPass.
-    However, this one is not well integrated into the BQSKit tool suite,
-    and as a result lacks customizability.
-
-    As we continue to port over the features that make this implementation
-    slightly higher quality, we will deprecate this pass.
-    """
-
-    def synthesize(self, utry: UnitaryMatrix, data: dict[str, Any]) -> Circuit:
-        # Set leap options
-        opts = options.Options()
-        opts.target = utry.numpy
-        opts.verbosity = 0
-        opts.write_to_stdout = False
-        opts.reoptimize_size = 7
-        opts.parallelizer = parallelizers.SequentialParallelizer
-
-        # Run the LEAP algorithm and optimizations
-        intermediate = leap_compiler.LeapCompiler().compile(opts)
-        pp = post_processing.LEAPReoptimizing_PostProcessor()
-        output = pp.post_process_circuit(intermediate, opts)
-
-        # Translate back to BQSKit through QASM
-        output = assemblers.ASSEMBLER_IBMOPENQASM.assemble(output)
-        return OPENQASM2Language().decode(output)
diff --git a/bqskit/passes/synthesis/qfast.py b/bqskit/passes/synthesis/qfast.py
index 8cbd658cf..e906a2304 100644
--- a/bqskit/passes/synthesis/qfast.py
+++ b/bqskit/passes/synthesis/qfast.py
@@ -66,7 +66,7 @@ def __init__(
 
             max_depth (int): The maximum number of gates to append without
                 success before termination. If left as None it will default
-                 to unlimited. (Default: None)
+                to unlimited. (Default: None)
 
         Raises:
             ValueError: If max_depth is nonpositive.
diff --git a/bqskit/passes/synthesis/qpredict.py b/bqskit/passes/synthesis/qpredict.py
index 0f7488080..fee8e0443 100644
--- a/bqskit/passes/synthesis/qpredict.py
+++ b/bqskit/passes/synthesis/qpredict.py
@@ -75,7 +75,7 @@ def __init__(
 
             max_depth (int): The maximum number of gates to append without
                 success before termination. If left as None it will default
-                 to unlimited. (Default: None)
+                to unlimited. (Default: None)
 
             instantiate_options (dict[str: Any]): Options passed directly
                 to circuit.instantiate when instantiating circuit
diff --git a/bqskit/passes/synthesis/qsearch.py b/bqskit/passes/synthesis/qsearch.py
index af1547166..dbf2338f9 100644
--- a/bqskit/passes/synthesis/qsearch.py
+++ b/bqskit/passes/synthesis/qsearch.py
@@ -4,8 +4,7 @@
 import logging
 from typing import Any
 
-from dask.distributed import as_completed
-from dask.distributed import worker_client
+from dask.distributed import get_client
 
 from bqskit.ir.circuit import Circuit
 from bqskit.ir.opt.cost.functions import HilbertSchmidtResidualsGenerator
@@ -41,6 +40,8 @@ def __init__(
         success_threshold: float = 1e-10,
         cost: CostFunctionGenerator = HilbertSchmidtResidualsGenerator(),
         max_layer: int | None = None,
+        store_partial_solutions: bool = False,
+        partials_per_depth: int = 25,
         instantiate_options: dict[str, Any] = {},
     ) -> None:
         """
@@ -65,7 +66,14 @@ def __init__(
 
             max_layer (int): The maximum number of layers to append without
                 success before termination. If left as None it will default
-                 to unlimited. (Default: None)
+                to unlimited. (Default: None)
+
+            store_partial_solutions (bool): Whether to store partial solutions
+                at different depths inside of the data dict. (Default: False)
+
+            partials_per_depth (int): The maximum number of partials
+                to store per search depth. No effect if
+                `store_partial_solutions` is False. (Default: 25)
 
             instantiate_options (dict[str: Any]): Options passed directly
                 to circuit.instantiate when instantiating circuit
@@ -121,6 +129,8 @@ def __init__(
         self.max_layer = max_layer
         self.instantiate_options: dict[str, Any] = {'cost_fn_gen': self.cost}
         self.instantiate_options.update(instantiate_options)
+        self.store_partial_solutions = store_partial_solutions
+        self.partials_per_depth = partials_per_depth
 
     def synthesize(self, utry: UnitaryMatrix, data: dict[str, Any]) -> Circuit:
         """Synthesize `utry`, see :class:`SynthesisPass` for more."""
@@ -140,6 +150,8 @@ def synthesize(self, utry: UnitaryMatrix, data: dict[str, Any]) -> Circuit:
         search_data['best_dist'] = self.cost.calc_cost(initial_layer, utry)
         search_data['best_circ'] = initial_layer
         search_data['best_layer'] = 0
+        if self.store_partial_solutions:
+            search_data['psols'] = {}
         _logger.info(
             'Search started, initial layer has cost: %e.' %
             search_data['best_dist'],
@@ -149,36 +161,50 @@ def synthesize(self, utry: UnitaryMatrix, data: dict[str, Any]) -> Circuit:
             _logger.info('Successful synthesis.')
             return initial_layer
 
-        while not frontier.empty():
-            top_circuit, layer = frontier.pop()
-
-            # Generate successors and evaluate each
-            successors = self.layer_gen.gen_successors(top_circuit, data)
-
-            if 'executor' in data:  # In Parallel
-                with worker_client() as client:
-                    futures = client.map(
-                        Circuit.instantiate,
-                        successors,
-                        pure=False,
-                        target=utry,
-                        **instantiate_options,
-                    )
-                    for _, circuit in as_completed(futures, with_results=True):
-                        if self.evaluate_node(
-                            circuit,
-                            utry,
-                            data,
-                            frontier,
-                            layer,
-                            search_data,
-                        ):
-                            for future in futures:
-                                if not future.done:
-                                    client.cancel(future)
-                            return circuit
-
-            else:  # Sequentially
+        if 'executor' in data:  # In Parallel
+            client = get_client()
+            while not frontier.empty():
+                top_circuit, layer = frontier.pop()
+
+                # Generate successors and evaluate each
+                successors = self.layer_gen.gen_successors(top_circuit, data)
+
+                # Submit instantiate jobs
+                futures = self.batched_instantiate(
+                    successors,
+                    utry,
+                    client,
+                    **self.instantiate_options,
+                )
+
+                # Wait for and gather results
+                circuits = self.gather_best_results(
+                    futures,
+                    client,
+                    self.cost.calc_cost,
+                    utry,
+                )
+
+                for circuit in circuits:
+                    if self.evaluate_node(
+                        circuit,
+                        utry,
+                        data,
+                        frontier,
+                        layer,
+                        search_data,
+                    ):
+                        if self.store_partial_solutions:
+                            data['psols'] = search_data['psols']
+                        return circuit
+
+        else:  # Sequentially
+            while not frontier.empty():
+                top_circuit, layer = frontier.pop()
+
+                # Generate successors and evaluate each
+                successors = self.layer_gen.gen_successors(top_circuit, data)
+
                 for circuit in successors:
                     circuit.instantiate(utry, **instantiate_options)
                     if self.evaluate_node(
@@ -189,6 +215,8 @@ def synthesize(self, utry: UnitaryMatrix, data: dict[str, Any]) -> Circuit:
                         layer,
                         search_data,
                     ):
+                        if self.store_partial_solutions:
+                            data['psols'] = search_data['psols']
                         return circuit
 
         _logger.info('Frontier emptied.')
@@ -199,6 +227,8 @@ def synthesize(self, utry: UnitaryMatrix, data: dict[str, Any]) -> Circuit:
                 else 's', search_data['best_dist'],
             ),
         )
+        if self.store_partial_solutions:
+            data['psols'] = search_data['psols']
 
         return search_data['best_circ']
 
@@ -226,6 +256,16 @@ def evaluate_node(
             search_data['best_circ'] = circuit
             search_data['best_layer'] = layer
 
+        if self.store_partial_solutions:
+            if layer not in search_data['psols']:
+                search_data['psols'][layer] = []
+
+            search_data['psols'][layer].append((circuit.copy(), dist))
+
+            if len(search_data['psols'][layer]) > self.partials_per_depth:
+                search_data['psols'][layer].sort(key=lambda x: x[1])
+                del search_data['psols'][layer][-1]
+
         if self.max_layer is None or layer + 1 < self.max_layer:
             frontier.add(circuit, layer + 1)
 
diff --git a/bqskit/passes/synthesis/synthesis.py b/bqskit/passes/synthesis/synthesis.py
index 9d78d1f9a..0b419bf9f 100644
--- a/bqskit/passes/synthesis/synthesis.py
+++ b/bqskit/passes/synthesis/synthesis.py
@@ -4,9 +4,18 @@
 import logging
 from abc import abstractmethod
 from typing import Any
+from typing import Callable
+
+from dask.distributed import as_completed
+from dask.distributed import Client
+from dask.distributed import Future
+from dask.distributed import rejoin
+from dask.distributed import secede
 
 from bqskit.compiler.basepass import BasePass
 from bqskit.ir.circuit import Circuit
+from bqskit.qis.state.state import StateLike
+from bqskit.qis.unitary.unitarymatrix import UnitaryLike
 from bqskit.qis.unitary.unitarymatrix import UnitaryMatrix
 
 _logger = logging.getLogger(__name__)
@@ -46,3 +55,154 @@ def run(self, circuit: Circuit, data: dict[str, Any] = {}) -> None:
 
         target_utry = circuit.get_unitary()
         circuit.become(self.synthesize(target_utry, data))
+
+    def batched_instantiate(
+        self,
+        circuits: list[Circuit],
+        target: UnitaryLike | StateLike,
+        client: Client,
+        **kwargs: Any,
+    ) -> list[list[Future]]:
+        """
+        Batch instantiate `circuits`.
+
+        Args:
+            circuits (list[Circuit]): The circuit batch to instantiate.
+
+            target (UnitaryLike | StateLike): The instantiation target.
+
+            client (Client): The Dask client used to submit jobs.
+
+            kwargs (Any): Other keyword arguments are passed directly to
+                instantiate calls.
+
+        Returns:
+            (list[list[Future]]): The Dask Futures corresponding to the
+                the instantiate jobs submitted. This returns a double list,
+                where the other list is indexed by input circuits, and the
+                inner list is indexed by multistart.
+        """
+        multistarts = 1
+        if 'multistarts' in kwargs:
+            multistarts = kwargs['multistarts']
+            kwargs['multistarts'] = 1
+
+        futures: list[list[Future]] = []
+
+        client.scatter(circuits)
+
+        for circuit in circuits:
+            futures.append([])
+            for i in range(multistarts):
+                futures[-1].append(
+                    client.submit(
+                        Circuit.instantiate,
+                        circuit,
+                        pure=False,
+                        target=target,
+                        **kwargs,
+                    ),
+                )
+
+        return futures
+
+    def gather_best_results(
+        self,
+        futures: list[list[Future]],
+        client: Client,
+        fn: Callable[..., float],
+        *args: Any,
+        **kwargs: Any,
+    ) -> list[Circuit]:
+        """
+        Gather best results from a `batched_instantiate` call.
+
+        Args:
+            futures (list[list[Future]]): The futures return from a
+                `batched_instantiate` call.
+
+            client (Client): The Dask client used to submit the jobs.
+
+            fn (Callable[..., float]): The function used to sort
+                instantiated circuits. This should take a circuit as
+                the first parameter.
+
+            args (Any): Arguments passed directly to fn.
+
+            kwargs (Any): Keyword arguments passed directly to fn.
+
+        Returns:
+            (list[Circuit]): The resulting circuits. There is one circuit
+                for each future list, i.e., `len(output) == len(futures)`.
+        """
+
+        score_futures: list[list[Future]] = []
+        for future_list in futures:
+            score_futures.append([])
+            for future in future_list:
+                score_futures[-1].append(
+                    client.submit(
+                        fn,
+                        future,
+                        *args,
+                        **kwargs,
+                    ),
+                )
+
+        flat_score_list = []
+        for future_list in score_futures:
+            flat_score_list.extend(future_list)
+
+        secede()
+        client.gather(flat_score_list)
+        rejoin()
+
+        best_circuit_futures = []
+        for i, score_list in enumerate(score_futures):
+            scores = client.gather(score_list)
+            best_index = scores.index(min(scores))
+            best_circuit_futures.append(futures[i][best_index])
+
+        return client.gather(best_circuit_futures)
+
+    def gather_first_results(
+        self,
+        futures: list[list[Future]],
+        client: Client,
+    ) -> list[Circuit]:
+        """
+        Gather first results from a `batched_instantiate` call.
+
+        Args:
+            futures (list[list[Future]]): The futures return from a
+                `batched_instantiate` call.
+
+            client (Client): The Dask client used to submit the jobs.
+
+        Returns:
+            (list[Circuit]): The resulting circuits. There is one circuit
+                for each future list, i.e., `len(output) == len(futures)`.
+                This call will return the first multistart instantiate
+                job for each circuit.
+        """
+        all_futures = []
+        future_to_index_map = {}
+        for index, future_list in enumerate(futures):
+            all_futures.extend(future_list)
+            for future in future_list:
+                future_to_index_map[future] = index
+
+        secede()
+        circuits = []
+        indices_seen = set()
+        for future, circuit in as_completed(all_futures, with_results=True):
+            index = future_to_index_map[future]
+            if index not in indices_seen:
+                indices_seen.add(index)
+                circuits.append((index, circuit))
+                for other_future in futures[index]:
+                    if not other_future.done():
+                        client.cancel(other_future)
+        rejoin()
+
+        return [c for _, c in sorted(circuits, key=lambda x: x[0])]
diff --git a/bqskit/passes/util/__init__.py b/bqskit/passes/util/__init__.py
index b5dbccfac..ea1776716 100644
--- a/bqskit/passes/util/__init__.py
+++ b/bqskit/passes/util/__init__.py
@@ -2,6 +2,7 @@
 from __future__ import annotations
 
 from bqskit.passes.util.compress import CompressPass
+from bqskit.passes.util.conversion import BlockConversionPass
 from bqskit.passes.util.converttou3 import ToU3Pass
 from bqskit.passes.util.random import SetRandomSeedPass
 from bqskit.passes.util.record import RecordStatsPass
@@ -15,4 +16,5 @@
     'UnfoldPass',
     'UpdateDataPass',
     'ToU3Pass',
+    'BlockConversionPass',
 ]
diff --git a/bqskit/passes/util/conversion.py b/bqskit/passes/util/conversion.py
new file mode 100644
index 000000000..60b75e2b2
--- /dev/null
+++ b/bqskit/passes/util/conversion.py
@@ -0,0 +1,114 @@
+"""This module implements the BlockConversionPass."""
+from __future__ import annotations
+
+import logging
+from typing import Any
+
+from bqskit.compiler.basepass import BasePass
+from bqskit.ir.circuit import Circuit
+from bqskit.ir.gates import CircuitGate
+from bqskit.ir.gates import ConstantUnitaryGate
+from bqskit.ir.gates import VariableUnitaryGate
+from bqskit.ir.point import CircuitPoint
+
+_logger = logging.getLogger(__name__)
+
+
+class BlockConversionPass(BasePass):
+    """
+    Converts blocks of one type to another type.
+
+    Blocks are either described by a constant or variable unitary gate or as a
+    circuit gate. Often during the flow of compilatin we will need to convert
+    them from one form to another for a future pass.
+    """
+
+    def __init__(
+        self,
+        convert_target: str,
+        convert_variable: bool = True,
+        convert_constant: bool = True,
+        convert_circuitgates: bool = True,
+    ):
+        """
+        Construct a BlockConversionPass.
+
+        Args:
+            convert_target (str): Either `variable` or `constant`.
+                Blocks will be converted to the form described here. If
+                this is `variable` all gates will be converted to
+                `VariableUnitaryGate` s. If this is `constant` blocks
+                will be converted to `ConstantUnitaryGate` s. Blocks
+                cannot be converted to circuit gates, that can be caried
+                out by synthesis.
+
+            convert_variable (bool): If this is true, this will replace
+                VariableUnitaryGate's in the circuit with one's specified
+                in convert_target.
+
+            convert_constant (bool): If this is true, this will replace
+                ConstantUnitaryGate's in the circuit with one's specified
+                in convert_target.
+
+            convert_circuitgates (bool): If this is true, this will replace
+                CircuitGate's in the circuit with one's specified
+                in convert_target. The subcircuit information captured
+                in the circuit gate will be lost.
+        """
+        self.convert_variable = convert_variable
+        self.convert_constant = convert_constant
+        self.convert_circuitgates = convert_circuitgates
+
+        if convert_target == 'variable':
+            self.convert_target = 'variable'
+            self.convert_variable = False
+        elif convert_target == 'constant':
+            self.convert_target = 'constant'
+            self.convert_constant = False
+        else:
+            raise ValueError('Unexpected input for conversion target.')
+
+    def run(self, circuit: Circuit, data: dict[str, Any] = {}) -> None:
+        """Perform the pass's operation, see :class:`BasePass` for more."""
+
+        # Variable -> Constant
+        if self.convert_variable and self.convert_target == 'constant':
+            _logger.debug('Converting variable gates to constant gates.')
+
+            for cycle, op in circuit.operations_with_cycles():
+                if isinstance(op.gate, VariableUnitaryGate):
+                    cgate = ConstantUnitaryGate(op.get_unitary(), op.radixes)
+                    point = CircuitPoint(cycle, op.location[0])
+                    circuit.replace_gate(point, cgate, op.location)
+
+        # CircuitGates -> Constant
+        if self.convert_circuitgates and self.convert_target == 'constant':
+            _logger.debug('Converting circuit gates to constant gates.')
+
+            for cycle, op in circuit.operations_with_cycles():
+                if isinstance(op.gate, CircuitGate):
+                    cgate = ConstantUnitaryGate(op.get_unitary(), op.radixes)
+                    point = CircuitPoint(cycle, op.location[0])
+                    circuit.replace_gate(point, cgate, op.location)
+
+        # Constant -> Variable
+        if self.convert_constant and self.convert_target == 'variable':
+            _logger.debug('Converting constant gates to variable gates.')
+
+            for cycle, op in circuit.operations_with_cycles():
+                if isinstance(op.gate, ConstantUnitaryGate):
+                    params = VariableUnitaryGate.get_params(op.get_unitary())
+                    vgate = VariableUnitaryGate(op.num_qudits, op.radixes)
+                    point = CircuitPoint(cycle, op.location[0])
+                    circuit.replace_gate(point, vgate, op.location, params)
+
+        # CircuitGates -> Variable
+        if self.convert_constant and self.convert_target == 'variable':
+            _logger.debug('Converting circuit gates to variable gates.')
+
+            for cycle, op in circuit.operations_with_cycles():
+                if isinstance(op.gate, CircuitGate):
+                    params = VariableUnitaryGate.get_params(op.get_unitary())
+                    vgate = VariableUnitaryGate(op.num_qudits, op.radixes)
+                    point = CircuitPoint(cycle, op.location[0])
+                    circuit.replace_gate(point, vgate, op.location, params)
diff --git a/bqskit/passes/util/record.py b/bqskit/passes/util/record.py
index d384852f3..d2574d052 100644
--- a/bqskit/passes/util/record.py
+++ b/bqskit/passes/util/record.py
@@ -1,11 +1,14 @@
 """This module implements the RecordStatsPass class."""
 from __future__ import annotations
 
+import logging
 from typing import Any
 
 from bqskit.compiler.basepass import BasePass
 from bqskit.ir.circuit import Circuit
 
+_logger = logging.getLogger(__name__)
+
 
 class RecordStatsPass(BasePass):
     """
@@ -32,3 +35,10 @@ def run(self, circuit: Circuit, data: dict[str, Any] = {}) -> None:
             data[self.key] = []
 
         data[self.key].append(stats)
+
+        _logger.info('Circuit Current Stats')
+        _logger.info(f'{circuit.num_cycles} Cycles')
+        _logger.info(f'{circuit.num_operations} Operation Count')
+        _logger.info(f'{circuit.coupling_graph} Connectivity')
+        for gate in circuit.gate_set:
+            _logger.info(f'{circuit.count(gate)} {gate.name} Count')
diff --git a/bqskit/qis/pauli.py b/bqskit/qis/pauli.py
index 28fe789c8..46becc42f 100644
--- a/bqskit/qis/pauli.py
+++ b/bqskit/qis/pauli.py
@@ -6,16 +6,20 @@
 from typing import Iterator
 from typing import overload
 from typing import Sequence
+from typing import TYPE_CHECKING
 
 import numpy as np
+import numpy.typing as npt
 
-from bqskit.qis.unitary.unitary import RealVector
 from bqskit.utils.typing import is_integer
 from bqskit.utils.typing import is_numeric
 from bqskit.utils.typing import is_sequence
 
+if TYPE_CHECKING:
+    from bqskit.qis.unitary.unitary import RealVector
 
-class PauliMatrices(Sequence[np.ndarray]):
+
+class PauliMatrices(Sequence[npt.NDArray[np.complex128]]):
     """
     Pauli group of matrices.
 
@@ -99,39 +103,44 @@ def __init__(self, num_qudits: int) -> None:
             for pauli_n_1, pauli_1 in matrices:
                 self.paulis.append(np.kron(pauli_n_1, pauli_1))
 
-    def __iter__(self) -> Iterator[np.ndarray]:
+    def __iter__(self) -> Iterator[npt.NDArray[np.complex128]]:
         return self.paulis.__iter__()
 
     @overload
-    def __getitem__(self, index: int) -> np.ndarray:
+    def __getitem__(self, index: int) -> npt.NDArray[np.complex128]:
         ...
 
     @overload
-    def __getitem__(self, index: slice) -> list[np.ndarray]:
+    def __getitem__(self, index: slice) -> list[npt.NDArray[np.complex128]]:
         ...
 
-    def __getitem__(self, index: int | slice) -> np.ndarray | list[np.ndarray]:
+    def __getitem__(
+        self,
+        index: int | slice,
+    ) -> npt.NDArray[np.complex128] | list[npt.NDArray[np.complex128]]:
         return self.paulis[index]
 
     def __len__(self) -> int:
         return len(self.paulis)
 
     @property
-    def numpy(self) -> np.ndarray:
+    def numpy(self) -> npt.NDArray[np.complex128]:
         """The NumPy array holding the pauli matrices."""
         return np.array(self.paulis)
 
     def __array__(
         self,
         dtype: np.typing.DTypeLike = np.complex128,
-    ) -> np.ndarray:
+    ) -> npt.NDArray[np.complex128]:
         """Implements NumPy API for the PauliMatrices class."""
         if dtype != np.complex128:
             raise ValueError('PauliMatrices only supports Complex128 dtype.')
 
         return np.array(self.paulis, dtype)
 
-    def get_projection_matrices(self, q_set: Iterable[int]) -> list[np.ndarray]:
+    def get_projection_matrices(
+            self, q_set: Iterable[int],
+    ) -> list[npt.NDArray[np.complex128]]:
         """
         Return the Pauli matrices that act only on qubits in `q_set`.
 
@@ -171,7 +180,7 @@ def get_projection_matrices(self, q_set: Iterable[int]) -> list[np.ndarray]:
 
         return pauli_n_qubit
 
-    def dot_product(self, alpha: RealVector) -> np.ndarray:
+    def dot_product(self, alpha: RealVector) -> npt.NDArray[np.complex128]:
         """
         Computes the standard dot product of `alpha` with the paulis.
 
@@ -200,7 +209,9 @@ def dot_product(self, alpha: RealVector) -> np.ndarray:
         return np.array(np.sum([a * s for a, s in zip(alpha, self.paulis)], 0))
 
     @staticmethod
-    def from_string(pauli_string: str) -> np.ndarray | list[np.ndarray]:
+    def from_string(
+            pauli_string: str,
+    ) -> npt.NDArray[np.complex128] | list[npt.NDArray[np.complex128]]:
         """
         Construct pauli matrices from a string description.
 
diff --git a/bqskit/qis/state/state.py b/bqskit/qis/state/state.py
index 255030822..b653628e3 100644
--- a/bqskit/qis/state/state.py
+++ b/bqskit/qis/state/state.py
@@ -8,6 +8,7 @@
 from typing import Union
 
 import numpy as np
+import numpy.typing as npt
 from numpy.lib.mixins import NDArrayOperatorsMixin
 from scipy.stats import unitary_group
 
@@ -95,7 +96,7 @@ def __init__(
         self._dim = dim
 
     @property
-    def numpy(self) -> np.ndarray:
+    def numpy(self) -> npt.NDArray[np.complex128]:
         """The NumPy array holding the vector."""
         return self._vec
 
@@ -132,7 +133,9 @@ def __iter__(self) -> Iterator[np.complex128]:
         """An iterator that iterates through the elements of the vector."""
         return self._vec.__iter__()
 
-    def __getitem__(self, index: Any) -> np.complex128 | np.ndarray:
+    def __getitem__(
+            self, index: Any,
+    ) -> np.complex128 | npt.NDArray[np.complex128]:
         """Implements NumPy API for the StateVector class."""
         return self._vec[index]
 
@@ -217,7 +220,7 @@ def __eq__(self, other: object) -> bool:
     def __array__(
         self,
         dtype: np.typing.DTypeLike = np.complex128,
-    ) -> np.ndarray:
+    ) -> npt.NDArray[np.complex128]:
         """Implements NumPy API for the StateVector class."""
         if dtype != np.complex128:
             raise ValueError('StateVector only supports Complex128 dtype.')
@@ -228,15 +231,15 @@ def __array_ufunc__(
         self,
         ufunc: np.ufunc,
         method: str,
-        *inputs: np.ndarray,
+        *inputs: npt.NDArray[Any],
         **kwargs: Any,
-    ) -> StateVector | np.ndarray:
+    ) -> StateVector | npt.NDArray[np.complex128]:
         """Implements NumPy API for the StateVector class."""
         if method != '__call__':
             return NotImplemented
 
         non_state_involved = False
-        args: list[np.ndarray] = []
+        args: list[npt.NDArray[Any]] = []
         for input in inputs:
             if isinstance(input, StateVector):
                 args.append(input.numpy)
diff --git a/bqskit/qis/unitary/differentiable.py b/bqskit/qis/unitary/differentiable.py
index 5453c6d6f..b2348d008 100644
--- a/bqskit/qis/unitary/differentiable.py
+++ b/bqskit/qis/unitary/differentiable.py
@@ -4,6 +4,7 @@
 import abc
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.qis.unitary.unitary import RealVector
 from bqskit.qis.unitary.unitary import Unitary
@@ -20,7 +21,7 @@ class DifferentiableUnitary(Unitary):
     """
 
     @abc.abstractmethod
-    def get_grad(self, params: RealVector = []) -> np.ndarray:
+    def get_grad(self, params: RealVector = []) -> npt.NDArray[np.complex128]:
         """
         Return the gradient for the unitary map as an np.ndarray.
 
@@ -44,7 +45,7 @@ def get_grad(self, params: RealVector = []) -> np.ndarray:
     def get_unitary_and_grad(
         self,
         params: RealVector = [],
-    ) -> tuple[UnitaryMatrix, np.ndarray]:
+    ) -> tuple[UnitaryMatrix, npt.NDArray[np.complex128]]:
         """
         Return a tuple combining the outputs of `get_unitary` and `get_grad`.
 
diff --git a/bqskit/qis/unitary/optimizable.py b/bqskit/qis/unitary/optimizable.py
index 15e3beeeb..b5dff6fd7 100644
--- a/bqskit/qis/unitary/optimizable.py
+++ b/bqskit/qis/unitary/optimizable.py
@@ -4,6 +4,7 @@
 import abc
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.qis.unitary.unitary import Unitary
 from bqskit.utils.typing import is_square_matrix
@@ -20,7 +21,7 @@ class LocallyOptimizableUnitary(Unitary):
     """
 
     @abc.abstractmethod
-    def optimize(self, env_matrix: np.ndarray) -> list[float]:
+    def optimize(self, env_matrix: npt.NDArray[np.complex128]) -> list[float]:
         """
         Return the optimal parameters with respect to an environment matrix.
 
@@ -34,14 +35,14 @@ def optimize(self, env_matrix: np.ndarray) -> list[float]:
             \\mathit{self.get\\_unitary(params))}))
 
         Args:
-            env_matrix (np.ndarray): Optimize with respect to this matrix.
-                Has the same dimensions as this unitary.
+            env_matrix (npt.NDArray[np.complex128]): Optimize with respect
+                to this matrix. Has the same dimensions as this unitary.
 
         Returns:
             list[float]: The parameters that optimizes this unitary.
         """
 
-    def check_env_matrix(self, env_matrix: np.ndarray) -> None:
+    def check_env_matrix(self, env_matrix: npt.NDArray[np.complex128]) -> None:
         """Check to ensure the `env_matrix` is validly shaped."""
         if not is_square_matrix(env_matrix):
             raise TypeError('Expected a square matrix.')
diff --git a/bqskit/qis/unitary/unitarybuilder.py b/bqskit/qis/unitary/unitarybuilder.py
index 202ca4b01..ada56d67a 100644
--- a/bqskit/qis/unitary/unitarybuilder.py
+++ b/bqskit/qis/unitary/unitarybuilder.py
@@ -6,6 +6,7 @@
 from typing import Sequence
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.ir.location import CircuitLocation
 from bqskit.ir.location import CircuitLocationLike
@@ -75,7 +76,7 @@ def __init__(self, num_qudits: int, radixes: Sequence[int] = []) -> None:
 
         self._num_params = 0
         self._dim = int(np.prod(self.radixes))
-        self.tensor = np.identity(self.dim)
+        self.tensor = np.identity(self.dim, dtype=np.complex128)
         self.tensor = self.tensor.reshape(self.radixes * 2)
 
     def get_unitary(self, params: RealVector = []) -> UnitaryMatrix:
@@ -249,7 +250,9 @@ def apply_left(
         inv_perm = list(np.argsort(perm))
         self.tensor = self.tensor.transpose(inv_perm)
 
-    def calc_env_matrix(self, location: Sequence[int]) -> np.ndarray:
+    def calc_env_matrix(
+            self, location: Sequence[int],
+    ) -> npt.NDArray[np.complex128]:
         """
         Calculates the environment matrix w.r.t. the specified location.
 
diff --git a/bqskit/qis/unitary/unitarymatrix.py b/bqskit/qis/unitary/unitarymatrix.py
index 2d28f4b52..f83522d9f 100644
--- a/bqskit/qis/unitary/unitarymatrix.py
+++ b/bqskit/qis/unitary/unitarymatrix.py
@@ -2,28 +2,31 @@
 from __future__ import annotations
 
 import logging
-import os
 from typing import Any
 from typing import Sequence
 from typing import Union
 
 import numpy as np
+import numpy.typing as npt
 import scipy as sp
-
-if 'READTHEDOCS' not in os.environ:
-    from numpy.lib.mixins import NDArrayOperatorsMixin
-else:
-    class NDArrayOperatorsMixin:  # type: ignore
-        pass
 from scipy.stats import unitary_group
 
 from bqskit.qis.state.state import StateLike
 from bqskit.qis.state.state import StateVector
 from bqskit.qis.state.statemap import StateVectorMap
-from bqskit.qis.unitary.unitary import RealVector, Unitary
+from bqskit.qis.unitary.unitary import RealVector
+from bqskit.qis.unitary.unitary import Unitary
+from bqskit.utils.docs import building_docs
 from bqskit.utils.typing import is_integer
 from bqskit.utils.typing import is_square_matrix
 from bqskit.utils.typing import is_valid_radixes
+
+if not building_docs():
+    from numpy.lib.mixins import NDArrayOperatorsMixin
+else:
+    class NDArrayOperatorsMixin:  # type: ignore
+        pass
+
 _logger = logging.getLogger(__name__)
 
 
@@ -70,8 +73,8 @@ def __init__(
         """
 
         # Stop any actual logic when building documentation
-        if 'READTHEDOCS' in os.environ:
-            self._utry = np.array([])
+        if building_docs():
+            self._utry: npt.NDArray[np.complex128] = np.array([])
             return
 
         # Copy constructor
@@ -119,7 +122,7 @@ def __init__(
     _num_params = 0
 
     @property
-    def numpy(self) -> np.ndarray:
+    def numpy(self) -> npt.NDArray[np.complex128]:
         """The NumPy array holding the unitary."""
         return self._utry
 
@@ -258,7 +261,7 @@ def identity(dim: int, radixes: Sequence[int] = []) -> UnitaryMatrix:
 
     @staticmethod
     def closest_to(
-        M: np.ndarray,
+        M: npt.NDArray[np.complex128],
         radixes: Sequence[int] = [],
     ) -> UnitaryMatrix:
         """
@@ -330,7 +333,10 @@ def random(num_qudits: int, radixes: Sequence[int] = []) -> UnitaryMatrix:
     def __eq__(self, other: object) -> bool:
         """Check if `self` is approximately equal to `other`."""
         if isinstance(other, Unitary):
-            return np.allclose(self, other.get_unitary())
+            other_unitary = other.get_unitary()
+            if self.shape != other_unitary.shape:
+                return False
+            return np.allclose(self, other_unitary)
 
         if isinstance(other, np.ndarray):
             return np.allclose(self, other)
@@ -341,7 +347,9 @@ def save(self, filename: str) -> None:
         """Save the unitary to a file."""
         np.savetxt(filename, self.numpy)
 
-    def __getitem__(self, index: Any) -> np.complex128 | np.ndarray:
+    def __getitem__(
+            self, index: Any,
+    ) -> np.complex128 | npt.NDArray[np.complex128]:
         """Implements NumPy API for the StateVector class."""
         return self._utry[index]
 
@@ -400,7 +408,7 @@ def is_unitary(U: np.typing.ArrayLike, tol: float = 1e-8) -> bool:
     def __array__(
         self,
         dtype: np.typing.DTypeLike = np.complex128,
-    ) -> np.ndarray:
+    ) -> npt.NDArray[np.complex128]:
         """Implements NumPy API for the UnitaryMatrix class."""
         if dtype != np.complex128:
             raise ValueError('UnitaryMatrix only supports Complex128 dtype.')
@@ -411,15 +419,15 @@ def __array_ufunc__(
         self,
         ufunc: np.ufunc,
         method: str,
-        *inputs: np.ndarray,
+        *inputs: npt.NDArray[Any],
         **kwargs: Any,
-    ) -> UnitaryMatrix | np.ndarray:
+    ) -> UnitaryMatrix | npt.NDArray[np.complex128]:
         """Implements NumPy API for the UnitaryMatrix class."""
         if method != '__call__':
             return NotImplemented
 
         non_unitary_involved = False
-        args: list[np.ndarray] = []
+        args: list[npt.NDArray[Any]] = []
         for input in inputs:
             if isinstance(input, UnitaryMatrix):
                 args.append(input.numpy)
@@ -465,5 +473,9 @@ def __repr__(self) -> str:
         """Return the repr representation of the unitary."""
         return repr(self._utry)
 
+    def __hash__(self) -> int:
+        """Return the hash of the unitary."""
+        return hash((self._utry[0][0], self._utry[-1][-1], self.shape))
+
 
 UnitaryLike = Union[UnitaryMatrix, np.ndarray, Sequence[Sequence[Any]]]
diff --git a/bqskit/utils/__init__.py b/bqskit/utils/__init__.py
index bbfa88c21..a20474f4d 100644
--- a/bqskit/utils/__init__.py
+++ b/bqskit/utils/__init__.py
@@ -21,6 +21,8 @@
 
     math
     typing
+    random
+    docs
     test.strategies
     test.types
 """
diff --git a/bqskit/utils/docs.py b/bqskit/utils/docs.py
new file mode 100644
index 000000000..834f706b2
--- /dev/null
+++ b/bqskit/utils/docs.py
@@ -0,0 +1,14 @@
+"""This module implements helper methods for generating docs."""
+from __future__ import annotations
+
+import os
+
+
+def building_docs() -> bool:
+    """Return true if currently building documentations."""
+
+    if 'READTHEDOCS' in os.environ:
+        return True
+    if '__SPHINX_BUILD__' in os.environ:
+        return True
+    return False
diff --git a/bqskit/utils/math.py b/bqskit/utils/math.py
index 0e7ebbe83..00f85d7af 100644
--- a/bqskit/utils/math.py
+++ b/bqskit/utils/math.py
@@ -1,14 +1,19 @@
 """This module implements numerical functions."""
 from __future__ import annotations
 
+from typing import Any
+
 import numpy as np
+import numpy.typing as npt
 import scipy as sp
 
 from bqskit.qis.pauli import PauliMatrices
 from bqskit.qis.unitary.unitary import RealVector
 
 
-def dexpmv(M: np.ndarray, dM: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
+def dexpmv(
+    M: npt.NDArray[np.complex128], dM: npt.NDArray[np.complex128],
+) -> tuple[npt.NDArray[np.complex128], npt.NDArray[np.complex128]]:
     """
     Compute the Matrix exponential F = e^M and its derivative dF.
 
@@ -72,7 +77,10 @@ def dexpmv(M: np.ndarray, dM: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
     return F, dF
 
 
-def softmax(x: np.ndarray, beta: int = 20) -> np.ndarray:
+def softmax(
+    x: npt.NDArray[np.float64],
+    beta: int = 20,
+) -> npt.NDArray[np.float64]:
     """
     Computes the softmax of vector x.
 
@@ -90,7 +98,7 @@ def softmax(x: np.ndarray, beta: int = 20) -> np.ndarray:
     return exps / np.sum(exps)
 
 
-def dot_product(alpha: RealVector, sigma: RealVector) -> np.ndarray:
+def dot_product(alpha: RealVector, sigma: RealVector) -> npt.NDArray[Any]:
     """
     Computes the standard dot product of `alpha` with `sigma`.
 
@@ -107,7 +115,9 @@ def dot_product(alpha: RealVector, sigma: RealVector) -> np.ndarray:
     return np.array(np.sum([a * s for a, s in zip(alpha, sigma)], 0))
 
 
-def unitary_log_no_i(U: np.ndarray) -> np.ndarray:
+def unitary_log_no_i(
+        U: npt.NDArray[np.complex128],
+) -> npt.NDArray[np.complex128]:
     """
     Solves for H in U = e^{iH}
 
@@ -130,7 +140,7 @@ def unitary_log_no_i(U: np.ndarray) -> np.ndarray:
     return 0.5 * H0 + 0.5 * H0.conj().T
 
 
-def pauli_expansion(H: np.ndarray) -> np.ndarray:
+def pauli_expansion(H: npt.NDArray[np.complex128]) -> npt.NDArray[np.float64]:
     """
     Computes a Pauli expansion of the hermitian matrix H.
 
diff --git a/bqskit/utils/random.py b/bqskit/utils/random.py
new file mode 100644
index 000000000..e7d415b21
--- /dev/null
+++ b/bqskit/utils/random.py
@@ -0,0 +1,15 @@
+"""This module implements helper methods for random generators."""
+from __future__ import annotations
+
+import ctypes
+
+import numpy as np
+
+
+def seed_random_sources(seed: int) -> None:
+    """Seed the various sources of randomness BQSKit uses."""
+    # set rand() seed, used by Ceres
+    libc = ctypes.CDLL(None)
+    libc.srand(seed)
+    # set numpy seed
+    np.random.seed(seed)
diff --git a/bqskit/utils/typing.py b/bqskit/utils/typing.py
index 810cb0583..d082b536d 100644
--- a/bqskit/utils/typing.py
+++ b/bqskit/utils/typing.py
@@ -10,6 +10,7 @@
 from typing import Sized
 
 import numpy as np
+import numpy.typing as npt
 from typing_extensions import TypeGuard
 
 
@@ -221,7 +222,7 @@ def is_square_matrix(M: np.typing.ArrayLike) -> bool:
     return True
 
 
-def is_hermitian(H: np.ndarray, tol: float = 1e-8) -> bool:
+def is_hermitian(H: npt.NDArray[np.complex128], tol: float = 1e-8) -> bool:
     """Return true if H is a hermitian matrix."""
 
     if not is_square_matrix(H):
@@ -239,7 +240,7 @@ def is_hermitian(H: np.ndarray, tol: float = 1e-8) -> bool:
     return True
 
 
-def is_skew_hermitian(H: np.ndarray, tol: float = 1e-8) -> bool:
+def is_skew_hermitian(H: npt.NDArray[np.complex128], tol: float = 1e-8) -> bool:
     """Return true if H is a skew hermitian matrix."""
 
     if not is_square_matrix(H):
diff --git a/docs/conf.py b/docs/conf.py
index fa116d879..48512a433 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -12,13 +12,14 @@
 
 import os
 import sys
-sys.path.insert(0, os.path.abspath('../bqskit'))
+sys.path.insert(0, os.path.abspath('..'))
 
+os.environ['__SPHINX_BUILD__'] = 'True'
 
 # -- Project information -----------------------------------------------------
 
 project = 'BQSKit'
-copyright = '2021, Lawrence Berkeley National Laboratory'
+copyright = '2022, Lawrence Berkeley National Laboratory'
 author = 'BQSKit Development Team'
 
 
@@ -94,9 +95,9 @@
     'bqskitrs',
     'threadpoolctl',
     'dask',
-    'qsearch',
     'scipy',
     'numpy',
+    'qiskit',
     'lark-parser',
     'distributed',
 ]
diff --git a/docs/index.rst b/docs/index.rst
index 2267951dc..75ecf527f 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -36,3 +36,4 @@ and algorithm framework that can run many algorithms efficiently.
    source/passes
    source/qis
    source/utils
+   source/exec
diff --git a/docs/source/exec.rst b/docs/source/exec.rst
new file mode 100644
index 000000000..285e14ed0
--- /dev/null
+++ b/docs/source/exec.rst
@@ -0,0 +1,4 @@
+.. automodule:: bqskit.exec
+   :no-members:
+   :no-inherited-members:
+   :no-special-members:
diff --git a/examples/optimize.py b/examples/optimize.py
index 2a9fe57bc..d79bafb53 100644
--- a/examples/optimize.py
+++ b/examples/optimize.py
@@ -11,5 +11,7 @@
 task = CompilationTask.optimize(circuit)
 
 # Spawn a compiler and compile the task
-with Compiler() as compiler:
-    optimized_circuit = compiler.compile(task)
+if __name__ == '__main__':
+    with Compiler() as compiler:
+        optimized_circuit = compiler.compile(task)
+        print(optimized_circuit)
diff --git a/examples/qfast.py b/examples/qfast.py
index a67cff982..27dfb17ef 100644
--- a/examples/qfast.py
+++ b/examples/qfast.py
@@ -1,8 +1,6 @@
 """This script is contains a simple use case of the QFAST synthesis method."""
 from __future__ import annotations
 
-import logging
-
 from bqskit.compiler import CompilationTask
 from bqskit.compiler import Compiler
 from bqskit.ir import Circuit
@@ -14,9 +12,6 @@
 from bqskit.qis import UnitaryMatrix
 
 if __name__ == '__main__':
-    # Enable logging
-    logging.getLogger('bqskit').setLevel(logging.DEBUG)
-
     # Let's create a random 3-qubit unitary to synthesize and add it to a
     # circuit.
     circuit = Circuit.from_unitary(UnitaryMatrix.random(3))
@@ -33,3 +28,4 @@
     # Finally let's create create the compiler and execute the CompilationTask.
     with Compiler() as compiler:
         compiled_circuit = compiler.compile(task)
+        print(compiled_circuit)
diff --git a/examples/quest.py b/examples/quest.py
new file mode 100644
index 000000000..8ef008043
--- /dev/null
+++ b/examples/quest.py
@@ -0,0 +1,44 @@
+"""This script synthesizes and executes a circuit using the Quest algorithm."""
+from __future__ import annotations
+
+import logging
+
+from bqskit.compiler import CompilationTask
+from bqskit.compiler import Compiler
+from bqskit.exec.runners.quest import QuestRunner
+from bqskit.exec.runners.sim import SimulationRunner
+from bqskit.ir import Circuit
+from bqskit.passes import ForEachBlockPass
+from bqskit.passes import LEAPSynthesisPass
+from bqskit.passes import QFASTDecompositionPass
+from bqskit.passes import ScanningGateRemovalPass
+from bqskit.passes import UnfoldPass
+from bqskit.qis import UnitaryMatrix
+
+if __name__ == '__main__':
+    # Enable logging
+    logging.getLogger('bqskit').setLevel(logging.INFO)
+
+    # Let's create a random 4-qubit unitary to synthesize and add it to a
+    # circuit.
+    circuit = Circuit.from_unitary(UnitaryMatrix.random(4))
+
+    # We will now define the CompilationTask we want to run.
+    task = CompilationTask(
+        circuit, [
+            QFASTDecompositionPass(),
+            ForEachBlockPass([LEAPSynthesisPass(), ScanningGateRemovalPass()]),
+            UnfoldPass(),
+        ],
+    )
+
+    # Finally let's create create the compiler and execute the CompilationTask.
+    with Compiler() as compiler:
+        compiled_circuit = compiler.compile(task)
+
+        # Execute the circuit with Quest
+        quest_runner = QuestRunner(SimulationRunner(), compiler=compiler)
+        # Use IBMQRunner instead of SimulationRunner to run on a machine.
+        results = quest_runner.run(compiled_circuit)
+        print('Output Probability Distribution')
+        print(results.probs)
diff --git a/examples/synthesis.py b/examples/synthesis.py
index 3004ec4ee..ee1094c39 100644
--- a/examples/synthesis.py
+++ b/examples/synthesis.py
@@ -21,5 +21,7 @@
 task = CompilationTask.synthesize(toffoli)
 
 # Spawn a compiler and compile the task
-with Compiler() as compiler:
-    synthesized_circuit = compiler.compile(task)
+if __name__ == '__main__':
+    with Compiler() as compiler:
+        synthesized_circuit = compiler.compile(task)
+        print(synthesized_circuit)
diff --git a/examples/tutorials/Partitioning and Optimization.ipynb b/examples/tutorials/Partitioning and Optimization.ipynb
index 865a54959..4c6f72743 100644
--- a/examples/tutorials/Partitioning and Optimization.ipynb
+++ b/examples/tutorials/Partitioning and Optimization.ipynb
@@ -205,12 +205,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from bqskit.passes import OptimizedLEAPPass\n",
+    "from bqskit.passes import LEAPSynthesisPass\n",
     "\n",
     "\n",
     "task = CompilationTask(circuit, [\n",
     "    QuickPartitioner(4),\n",
-    "    ForEachBlockPass([OptimizedLEAPPass]),\n",
+    "    ForEachBlockPass([LEAPSynthesisPass()]),\n",
     "    UnfoldPass(),\n",
     "])\n",
     "\n",
diff --git a/setup.cfg b/setup.cfg
index 8998756ba..338e7f81d 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,6 +1,6 @@
 [metadata]
 name = bqskit
-version = 0.3.0
+version = 0.4.0
 description = Berkeley Quantum Synthesis Toolkit
 long_description = file: README.md
 long_description_content_type = text/markdown
@@ -19,6 +19,7 @@ classifiers =
     Programming Language :: Python :: 3.7
     Programming Language :: Python :: 3.8
     Programming Language :: Python :: 3.9
+    Programming Language :: Python :: 3.10
     Topic :: Scientific/Engineering
     Topic :: Scientific/Engineering :: Mathematics
     Topic :: Scientific/Engineering :: Physics
@@ -36,10 +37,9 @@ install_requires =
     dask[distributed]
     lark-parser
     numpy>=1.16.6
-    qsearch
     scipy>=1.4.1
     threadpoolctl
-    typing_extensions
+    typing-extensions
 python_requires = >=3.7, <4
 
 [options.packages.find]
@@ -69,7 +69,7 @@ pretty = true
 
 [tox:tox]
 minversion = 3.3.0
-envlist = pre-commit,py37,py38,py39
+envlist = pre-commit,py37,py38,py39,py310
 skip_missing_interpreters = true
 
 [testenv]


Test patch:
diff --git a/bqskit/utils/test/strategies.py b/bqskit/utils/test/strategies.py
index fd676be09..327ab174a 100644
--- a/bqskit/utils/test/strategies.py
+++ b/bqskit/utils/test/strategies.py
@@ -352,7 +352,7 @@ def circuit_locations(
     min_size: int = 1,
     max_size: int = 5,
 ) -> CircuitLocation:
-    """Hypothesis strategy for generating `CircuitLocation`s."""
+    """Hypothesis strategy for generating `CircuitLocation` s."""
     idxs = integers(min_qudit, max_qudit)
     locations = lists(idxs, min_size=min_size, max_size=max_size, unique=True)
     return CircuitLocation(draw(locations))
@@ -366,7 +366,7 @@ def circuit_location_likes(
     min_size: int = 1,
     max_size: int = 5,
 ) -> CircuitLocationLike:
-    """Hypothesis strategy for generating `CircuitLocationLike`s."""
+    """Hypothesis strategy for generating `CircuitLocationLike` s."""
     idxs = integers(min_qudit, max_qudit)
     locations = lists(idxs, min_size=min_size, max_size=max_size, unique=True)
     location = draw(locations)
diff --git a/bqskit/utils/test/types.py b/bqskit/utils/test/types.py
index 33a9396a5..f884c596b 100644
--- a/bqskit/utils/test/types.py
+++ b/bqskit/utils/test/types.py
@@ -472,7 +472,7 @@ def invalid_type_test(
     valids = []
     invalids = []
     for id, param in inspect.signature(func_to_test).parameters.items():
-        if param.annotation == inspect._empty:  # type: ignore
+        if param.annotation == inspect._empty:
             raise ValueError(
                 'Need type annotation to generate invalid type tests.',
             )
@@ -544,7 +544,7 @@ def valid_type_test(
 
     strategies = []
     for id, param in inspect.signature(func_to_test).parameters.items():
-        if param.annotation == inspect._empty:  # type: ignore
+        if param.annotation == inspect._empty:
             raise ValueError(
                 'Need type annotation to generate invalid type tests.',
             )
diff --git a/tests/conftest.py b/tests/conftest.py
index eb7741989..b3590945f 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -12,6 +12,7 @@
 from typing import Sequence
 
 import numpy as np
+import numpy.typing as npt
 import pytest
 from hypothesis import HealthCheck
 from hypothesis import settings
@@ -64,7 +65,7 @@
 
 settings.register_profile(
     'default',
-    suppress_health_check=(HealthCheck.too_slow,),
+    suppress_health_check=(HealthCheck.too_slow, HealthCheck.filter_too_much),
 )
 settings.load_profile('default')
 
@@ -111,20 +112,16 @@
     'float_9': np.longdouble(1234.0),
     'float_10': np.float32(1234.0),
     'float_11': np.float64(1234.0),
-    # Needs typeshed sync
-    'complex_1': complex(0.0j),  # type: ignore
-    # Needs typeshed sync
-    'complex_2': complex(0.0 + 0.0j),    # type: ignore
+    'complex_1': complex(0.0j),
+    'complex_2': complex(0.0 + 0.0j),
     'complex_3': 1.0j,
     'complex_4': 1.0 + 1.0j,
     'complex_5': 1.0 - 1.0j,
     'complex_7': np.csingle(1234 + 1234j),
     'complex_8': np.cdouble(1234 + 1234j),
     'complex_9': np.clongdouble(1234 + 1234j),
-    # needs newer numpy
-    'complex_10': np.complex64(1234 + 1234j),  # type: ignore
-    # needs newer numpy
-    'complex_11': np.complex128(1234 + 1234j),  # type: ignore
+    'complex_10': np.complex64(1234 + 1234j),
+    'complex_11': np.complex128(1234 + 1234j),
     'bool_1': False,
     'bool_2': True,
     'bool_3': np.bool_(False),
@@ -148,7 +145,7 @@
         0.0j,
         1.23j,
         1.1 + 1.0j,
-        np.complex64(1.0 + 1.0j),  # type: ignore # needs newer numpy
+        np.complex64(1.0 + 1.0j),
     ],
     'seq-complex_4': [1.234e12j] * 10,
     'seq-bool_1': [],
@@ -455,7 +452,7 @@ def toffoli_unitary() -> UnitaryMatrix:
 
 
 @pytest.fixture
-def toffoli_unitary_np() -> np.ndarray:
+def toffoli_unitary_np() -> npt.NDArray[np.complex128]:
     return TOFFOLI
 
 
@@ -465,11 +462,11 @@ def swap_unitary() -> UnitaryMatrix:
 
 
 @pytest.fixture
-def swap_unitary_np() -> np.ndarray:
+def swap_unitary_np() -> npt.NDArray[np.complex128]:
     return SWAP
 
 
-def random_utry_gen(dims: int | Sequence[int]) -> np.ndarray:
+def random_utry_gen(dims: int | Sequence[int]) -> npt.NDArray[np.complex128]:
     """
     Generate a random unitary.
 
@@ -530,7 +527,8 @@ def invalid_utry_gen(dims: int | Sequence[int]) -> np.ndarray:
 
 
 @pytest.fixture
-def gen_invalid_utry_np() -> Callable[[int | Sequence[int]], np.ndarray]:
+def gen_invalid_utry_np(
+) -> Callable[[int | Sequence[int]], npt.NDArray[np.complex128]]:
     """Provide a method to generate random invalid unitaries."""
     return invalid_utry_gen
 
@@ -782,7 +780,7 @@ def circuit_gen(
 
         # 2a. Shuffle gates
         shuffled_gates = gateset
-        np.random.shuffle(shuffled_gates)
+        np.random.shuffle(shuffled_gates)  # type: ignore
 
         # 2b. Find first gate that is compatible
         for gate in shuffled_gates:
diff --git a/tests/exec/test_exec.py b/tests/exec/test_exec.py
new file mode 100644
index 000000000..777923f65
--- /dev/null
+++ b/tests/exec/test_exec.py
@@ -0,0 +1,18 @@
+from __future__ import annotations
+
+from bqskit.exec.runners.sim import SimulationRunner
+from bqskit.ir.circuit import Circuit
+from bqskit.ir.gates import CNOTGate
+from bqskit.ir.gates import HGate
+
+
+def test_sim() -> None:
+    circuit = Circuit(2)
+    circuit.append_gate(HGate(), 0)
+    circuit.append_gate(CNOTGate(), (0, 1))
+    results = SimulationRunner().run(circuit)
+    counts = results.get_counts(1024)
+    assert counts[0] == 512
+    assert counts[1] == 0
+    assert counts[2] == 0
+    assert counts[3] == 512
diff --git a/tests/ir/gates/parameterized/test_unitary.py b/tests/ir/gates/parameterized/test_unitary.py
index 1701db8ae..82a314b9b 100644
--- a/tests/ir/gates/parameterized/test_unitary.py
+++ b/tests/ir/gates/parameterized/test_unitary.py
@@ -51,3 +51,10 @@ def test_optimize(utry: UnitaryMatrix) -> None:
     vug = VariableUnitaryGate(utry.num_qudits, utry.radixes)
     params = vug.optimize(np.array(utry))
     assert vug.get_unitary(params).get_distance_from(utry.conj().T) < 1e-7
+
+
+@given(unitaries(2, (2,), 2))
+def test_is(utry: UnitaryMatrix) -> None:
+    vug1 = VariableUnitaryGate(utry.num_qudits, utry.radixes)
+    vug2 = VariableUnitaryGate(utry.num_qudits)
+    assert vug1 == vug2
diff --git a/tests/ir/gates/test_native.py b/tests/ir/gates/test_native.py
index a286c8d4d..07fce6072 100644
--- a/tests/ir/gates/test_native.py
+++ b/tests/ir/gates/test_native.py
@@ -9,6 +9,9 @@
 from bqskit.ir import Circuit as Circ
 from bqskit.ir.gate import Gate
 from bqskit.ir.gates import CNOTGate
+from bqskit.ir.gates import CRXGate
+from bqskit.ir.gates import CRYGate
+from bqskit.ir.gates import CRZGate
 from bqskit.ir.gates import RXGate
 from bqskit.ir.gates import RXXGate
 from bqskit.ir.gates import RYGate
@@ -23,6 +26,9 @@
 
 
 NATIVE_GATES: list[Gate] = [
+    CRXGate(),
+    CRYGate(),
+    CRZGate(),
     RXGate(),
     RYGate(),
     RZGate(),
diff --git a/tests/ir/test_interval.py b/tests/ir/test_interval.py
index 8ea841186..00bf49daf 100644
--- a/tests/ir/test_interval.py
+++ b/tests/ir/test_interval.py
@@ -99,17 +99,6 @@ def test_indices(interval: CycleInterval) -> None:
     assert is_sequence_of_int(indices)
 
 
-@given(cycle_intervals(), everything_except(int))
-def test_contains(interval: CycleInterval, other_type: Any) -> None:
-    for idx in interval:
-        assert idx in interval
-    assert interval.lower in interval
-    assert interval.upper in interval
-    assert interval.upper + 1 not in interval
-    assert interval.lower - 1 not in interval
-    assert other_type not in interval
-
-
 class TestCycleIntervalOverlaps:
 
     @valid_type_test(CycleInterval(0, 0).overlaps)
diff --git a/tests/passes/partitioning/test_scan.py b/tests/passes/partitioning/test_scan.py
index 476452348..79e474403 100644
--- a/tests/passes/partitioning/test_scan.py
+++ b/tests/passes/partitioning/test_scan.py
@@ -1,5 +1,7 @@
 from __future__ import annotations
 
+import pytest
+
 from bqskit.ir.circuit import Circuit
 from bqskit.ir.gates import CircuitGate
 from bqskit.ir.gates import CNOTGate
@@ -60,118 +62,6 @@ def test_run_r6(self, r6_qudit_circuit: Circuit) -> None:
         for cycle_index in range(r6_qudit_circuit.num_cycles):
             assert not r6_qudit_circuit._is_cycle_idle(cycle_index)
 
-    def test_corner_case_1(self) -> None:
-        circuit = Circuit(6)
-
-        circuit.append_gate(
-            ConstantUnitaryGate(
-                UnitaryMatrix.random(2),
-            ), [
-                1, 4,
-            ],
-        )
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [2])
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [3])
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [5])
-
-        circuit.append_gate(
-            ConstantUnitaryGate(
-                UnitaryMatrix.random(2),
-            ), [
-                5, 0,
-            ],
-        )
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [1])
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [3])
-
-        circuit.append_gate(
-            ConstantUnitaryGate(
-                UnitaryMatrix.random(4),
-            ), [
-                2, 0, 3, 5,
-            ],
-        )
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [1])
-
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [1])
-
-        utry = circuit.get_unitary()
-        ScanPartitioner(3).run(circuit, {})
-
-        assert all(
-            isinstance(op.gate, CircuitGate)
-            or len(op.location) > 3
-            for op in circuit
-        )
-        assert all(
-            not isinstance(op.gate, TaggedGate)
-            or not op.gate.tag != '__fold_placeholder__'
-            for op in circuit
-        )
-        assert circuit.get_unitary() == utry
-        for cycle_index in range(circuit.num_cycles):
-            assert not circuit._is_cycle_idle(cycle_index)
-
-    def test_corner_case_2(self) -> None:
-        circuit = Circuit(6)
-
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [0])
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [1])
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [5])
-
-        circuit.append_gate(
-            ConstantUnitaryGate(
-                UnitaryMatrix.random(2),
-            ), [
-                3, 0,
-            ],
-        )
-
-        circuit.append_gate(
-            ConstantUnitaryGate(
-                UnitaryMatrix.random(2),
-            ), [
-                5, 0,
-            ],
-        )
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [3])
-
-        circuit.append_gate(
-            ConstantUnitaryGate(
-                UnitaryMatrix.random(4),
-            ), [
-                4, 0, 1, 2,
-            ],
-        )
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [5])
-
-        circuit.append_gate(
-            ConstantUnitaryGate(
-                UnitaryMatrix.random(3),
-            ), [
-                5, 0, 1,
-            ],
-        )
-
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [1])
-
-        utry = circuit.get_unitary()
-        ScanPartitioner(3).run(circuit, {})
-
-        assert all(
-            isinstance(op.gate, CircuitGate)
-            or len(op.location) > 3
-            for op in circuit
-        )
-        assert all(
-            not isinstance(op.gate, TaggedGate)
-            or not op.gate.tag != '__fold_placeholder__'
-            for op in circuit
-        )
-        assert circuit.get_unitary() == utry
-        for cycle_index in range(circuit.num_cycles):
-            assert not circuit._is_cycle_idle(cycle_index)
-
     def test_corner_case_3(self) -> None:
         circuit = Circuit(6)
         circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [0])
@@ -287,138 +177,6 @@ def test_corner_case_4(self) -> None:
         for cycle_index in range(circuit.num_cycles):
             assert not circuit._is_cycle_idle(cycle_index)
 
-    def test_corner_case_5(self) -> None:
-        circuit = Circuit(6)
-
-        circuit.append_gate(
-            ConstantUnitaryGate(
-                UnitaryMatrix.random(2),
-            ), [
-                4, 1,
-            ],
-        )
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [2])
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [5])
-
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [1])
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [4])
-
-        circuit.append_gate(
-            ConstantUnitaryGate(
-                UnitaryMatrix.random(4),
-            ), [
-                4, 1, 2, 3,
-            ],
-        )
-
-        circuit.append_gate(
-            ConstantUnitaryGate(
-                UnitaryMatrix.random(2),
-            ), [
-                4, 1,
-            ],
-        )
-
-        circuit.append_gate(
-            ConstantUnitaryGate(
-                UnitaryMatrix.random(2),
-            ), [
-                5, 1,
-            ],
-        )
-
-        circuit.append_gate(
-            ConstantUnitaryGate(
-                UnitaryMatrix.random(2),
-            ), [
-                2, 1,
-            ],
-        )
-
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [2])
-
-        utry = circuit.get_unitary()
-        ScanPartitioner(3).run(circuit, {})
-
-        assert all(
-            isinstance(op.gate, CircuitGate)
-            or len(op.location) > 3
-            for op in circuit
-        )
-        assert all(
-            not isinstance(op.gate, TaggedGate)
-            or not op.gate.tag != '__fold_placeholder__'
-            for op in circuit
-        )
-        assert circuit.get_unitary() == utry
-        for cycle_index in range(circuit.num_cycles):
-            assert not circuit._is_cycle_idle(cycle_index)
-
-    def test_corner_case_6(self) -> None:
-        circuit = Circuit(6)
-
-        circuit.append_gate(
-            ConstantUnitaryGate(
-                UnitaryMatrix.random(2),
-            ), [
-                1, 5,
-            ],
-        )
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [2])
-
-        circuit.append_gate(
-            ConstantUnitaryGate(
-                UnitaryMatrix.random(3),
-            ), [
-                4, 0, 2,
-            ],
-        )
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [5])
-
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [0])
-        circuit.append_gate(
-            ConstantUnitaryGate(
-                UnitaryMatrix.random(2),
-            ), [
-                5, 1,
-            ],
-        )
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [2])
-
-        circuit.append_gate(
-            ConstantUnitaryGate(
-                UnitaryMatrix.random(4),
-            ), [
-                0, 2, 3, 4,
-            ],
-        )
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [1])
-
-        circuit.append_gate(
-            ConstantUnitaryGate(
-                UnitaryMatrix.random(2),
-            ), [
-                1, 5,
-            ],
-        )
-
-        utry = circuit.get_unitary()
-        ScanPartitioner(3).run(circuit, {})
-
-        assert all(
-            isinstance(op.gate, CircuitGate)
-            or len(op.location) > 3
-            for op in circuit
-        )
-        assert all(
-            not isinstance(op.gate, TaggedGate)
-            or not op.gate.tag != '__fold_placeholder__'
-            for op in circuit
-        )
-        assert circuit.get_unitary() == utry
-        for cycle_index in range(circuit.num_cycles):
-            assert not circuit._is_cycle_idle(cycle_index)
-
     def test_corner_case_7(self) -> None:
         circuit = Circuit(6)
 
@@ -558,97 +316,16 @@ def test_corner_case_8(self) -> None:
         for cycle_index in range(circuit.num_cycles):
             assert not circuit._is_cycle_idle(cycle_index)
 
-    def test_corner_case_9(self) -> None:
+    def test_bigger_block(self) -> None:
         circuit = Circuit(6)
 
         circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [0])
-        circuit.append_gate(
-            ConstantUnitaryGate(
-                UnitaryMatrix.random(3),
-            ), [
-                2, 1, 4,
-            ],
-        )
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [3])
-
-        circuit.append_gate(
-            ConstantUnitaryGate(
-                UnitaryMatrix.random(2),
-            ), [
-                5, 1,
-            ],
-        )
+        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [1])
         circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [2])
         circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [3])
         circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [4])
-
-        circuit.append_gate(
-            ConstantUnitaryGate(
-                UnitaryMatrix.random(4),
-            ), [
-                5, 1, 2, 4,
-            ],
-        )
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [3])
-
-        circuit.append_gate(
-            ConstantUnitaryGate(
-                UnitaryMatrix.random(2),
-            ), [
-                4, 1,
-            ],
-        )
-
-        utry = circuit.get_unitary()
-        ScanPartitioner(3).run(circuit, {})
-
-        assert all(
-            isinstance(op.gate, CircuitGate)
-            or len(op.location) > 3
-            for op in circuit
-        )
-        assert all(
-            not isinstance(op.gate, TaggedGate)
-            or not op.gate.tag != '__fold_placeholder__'
-            for op in circuit
-        )
-        assert circuit.get_unitary() == utry
-        for cycle_index in range(circuit.num_cycles):
-            assert not circuit._is_cycle_idle(cycle_index)
-
-    def test_corner_case_10(self) -> None:
-        circuit = Circuit(6)
-
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [0])
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [1])
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [4])
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [5])
-
-        circuit.append_gate(
-            ConstantUnitaryGate(
-                UnitaryMatrix.random(2),
-            ), [
-                1, 2,
-            ],
-        )
         circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [5])
 
-        circuit.append_gate(
-            ConstantUnitaryGate(
-                UnitaryMatrix.random(2),
-            ), [
-                4, 1,
-            ],
-        )
-
-        circuit.append_gate(
-            ConstantUnitaryGate(
-                UnitaryMatrix.random(2),
-            ), [
-                3, 1,
-            ],
-        )
-
         circuit.append_gate(
             ConstantUnitaryGate(
                 UnitaryMatrix.random(4),
@@ -656,22 +333,5 @@ def test_corner_case_10(self) -> None:
                 1, 2, 3, 4,
             ],
         )
-
-        circuit.append_gate(ConstantUnitaryGate(UnitaryMatrix.random(1)), [3])
-
-        utry = circuit.get_unitary()
-        ScanPartitioner(3).run(circuit, {})
-
-        assert all(
-            isinstance(op.gate, CircuitGate)
-            or len(op.location) > 3
-            for op in circuit
-        )
-        assert all(
-            not isinstance(op.gate, TaggedGate)
-            or not op.gate.tag != '__fold_placeholder__'
-            for op in circuit
-        )
-        assert circuit.get_unitary() == utry
-        for cycle_index in range(circuit.num_cycles):
-            assert not circuit._is_cycle_idle(cycle_index)
+        with pytest.raises(RuntimeError):
+            ScanPartitioner(3).run(circuit)
diff --git a/tests/passes/synthesis/test_qsearch.py b/tests/passes/synthesis/test_qsearch.py
new file mode 100644
index 000000000..5f96a1297
--- /dev/null
+++ b/tests/passes/synthesis/test_qsearch.py
@@ -0,0 +1,28 @@
+from __future__ import annotations
+
+from bqskit.ir.circuit import Circuit
+from bqskit.passes import QSearchSynthesisPass
+from bqskit.qis import UnitaryMatrix
+# from bqskit.ir.gates import CPIGate
+# from bqskit.ir.gates import U8Gate
+# from bqskit.passes import SimpleLayerGenerator
+
+
+class TestQSearch:
+
+    def test_small_qubit(self) -> None:
+        utry = UnitaryMatrix.random(2)
+        circuit = Circuit.from_unitary(utry)
+        qsearch = QSearchSynthesisPass()
+        qsearch.run(circuit)
+        dist = circuit.get_unitary().get_distance_from(utry)
+        assert dist <= 1e-5
+
+    # def test_small_qutrit(self) -> None:
+    #     utry = UnitaryMatrix.random(2, [3, 3])
+    #     circuit = Circuit.from_unitary(utry)
+    #     layer_gen = SimpleLayerGenerator(CPIGate(), U8Gate())
+    #     qsearch = QSearchSynthesisPass(layer_generator=layer_gen)
+    #     qsearch.run(circuit)
+    #     dist = circuit.get_unitary().get_distance_from(utry)
+    #     assert dist <= qsearch.success_threshold
diff --git a/tests/qis/test_pauli.py b/tests/qis/test_pauli.py
index fee0df405..f5f8a8c8b 100644
--- a/tests/qis/test_pauli.py
+++ b/tests/qis/test_pauli.py
@@ -4,6 +4,7 @@
 from typing import Any
 
 import numpy as np
+import numpy.typing as npt
 import pytest
 from hypothesis import given
 from hypothesis.strategies import integers
@@ -528,7 +529,7 @@ def test_invalid_value(self, invalid_alpha: RealVector) -> None:
             ),
         ],
     )
-    def test_size_1(self, alpha: RealVector, prod: np.ndarray) -> None:
+    def test_size_1(self, alpha: RealVector, prod: npt.NDArray[Any]) -> None:
         assert np.allclose(PauliMatrices(1).dot_product(alpha), prod)
 
     @pytest.mark.parametrize(
@@ -609,7 +610,10 @@ def test_size_1(self, alpha: RealVector, prod: np.ndarray) -> None:
             ),
         ],
     )
-    def test_size_2(self, alpha: RealVector, prod: np.ndarray) -> None:
+    def test_size_2(
+        self, alpha: RealVector,
+        prod: npt.NDArray[np.complex128],
+    ) -> None:
         assert np.allclose(PauliMatrices(2).dot_product(alpha), prod)
 
 
@@ -681,7 +685,10 @@ def test_invalid_value(self, invalid_str: str) -> None:
             ('\t XY  ,,\n\r\t\t', np.kron(PauliMatrices.X, PauliMatrices.Y)),
         ],
     )
-    def test_single(self, pauli_str: str, pauli_mat: np.ndarray) -> None:
+    def test_single(
+        self, pauli_str: str,
+        pauli_mat: npt.NDArray[np.complex128],
+    ) -> None:
         assert isinstance(PauliMatrices.from_string(pauli_str), np.ndarray)
         assert np.allclose(
             np.array(PauliMatrices.from_string(pauli_str)),
@@ -788,7 +795,10 @@ def test_single(self, pauli_str: str, pauli_mat: np.ndarray) -> None:
             ),
         ],
     )
-    def test_multi(self, pauli_str: str, pauli_mats: np.ndarray) -> None:
+    def test_multi(
+        self, pauli_str: str,
+        pauli_mats: list[npt.NDArray[np.complex128]],
+    ) -> None:
         paulis = PauliMatrices.from_string(pauli_str)
         assert isinstance(paulis, list)
         assert all(isinstance(pauli, np.ndarray) for pauli in paulis)
diff --git a/tests/qis/unitary/test_meta.py b/tests/qis/unitary/test_meta.py
index d9d962b0f..d485dd65e 100644
--- a/tests/qis/unitary/test_meta.py
+++ b/tests/qis/unitary/test_meta.py
@@ -2,6 +2,7 @@
 from __future__ import annotations
 
 import numpy as np
+import numpy.typing as npt
 
 from bqskit.qis.unitary.differentiable import DifferentiableUnitary
 from bqskit.qis.unitary.meta import UnitaryMeta
@@ -17,7 +18,9 @@ class test_class(LocallyOptimizableUnitary):
             def get_unitary(self, p: RealVector = []) -> UnitaryMatrix:
                 pass
 
-            def optimize(self, env_matrix: np.ndarray) -> list[float]:
+            def optimize(
+                    self, env_matrix: npt.NDArray[np.complex128],
+            ) -> list[float]:
                 pass
 
         assert isinstance(test_class, UnitaryMeta)
@@ -28,7 +31,9 @@ class test_class(LocallyOptimizableUnitary):
             def get_unitary(self, p: RealVector = []) -> UnitaryMatrix:
                 pass
 
-            def optimize(self, env_matrix: np.ndarray) -> list[float]:
+            def optimize(
+                    self, env_matrix: npt.NDArray[np.complex128],
+            ) -> list[float]:
                 pass
 
             def is_locally_optimizable(self) -> bool:
@@ -42,7 +47,9 @@ class test_class(LocallyOptimizableUnitary):
             def get_unitary(self, p: RealVector = []) -> UnitaryMatrix:
                 pass
 
-            def optimize(self, env_matrix: np.ndarray) -> list[float]:
+            def optimize(
+                    self, env_matrix: npt.NDArray[np.complex128],
+            ) -> list[float]:
                 pass
 
             def is_locally_optimizable(self) -> bool:
@@ -59,7 +66,9 @@ class test_class(DifferentiableUnitary):
             def get_unitary(self, p: RealVector = []) -> UnitaryMatrix:
                 pass
 
-            def get_grad(self, params: RealVector = []) -> np.ndarray:
+            def get_grad(
+                self, params: RealVector = [],
+            ) -> npt.NDArray[np.complex128]:
                 pass
 
         assert isinstance(test_class, UnitaryMeta)
@@ -70,7 +79,9 @@ class test_class(DifferentiableUnitary):
             def get_unitary(self, p: RealVector = []) -> UnitaryMatrix:
                 pass
 
-            def get_grad(self, params: RealVector = []) -> np.ndarray:
+            def get_grad(
+                self, params: RealVector = [],
+            ) -> npt.NDArray[np.complex128]:
                 pass
 
             def is_differentiable(self) -> bool:
@@ -84,7 +95,9 @@ class test_class(DifferentiableUnitary):
             def get_unitary(self, p: RealVector = []) -> UnitaryMatrix:
                 pass
 
-            def get_grad(self, params: RealVector = []) -> np.ndarray:
+            def get_grad(
+                self, params: RealVector = [],
+            ) -> npt.NDArray[np.complex128]:
                 pass
 
             def is_differentiable(self) -> bool:
diff --git a/tests/qis/unitary/test_unitarymatrix.py b/tests/qis/unitary/test_unitarymatrix.py
index bfb9ceef7..89d4786a4 100644
--- a/tests/qis/unitary/test_unitarymatrix.py
+++ b/tests/qis/unitary/test_unitarymatrix.py
@@ -2,6 +2,7 @@
 from __future__ import annotations
 
 import numpy as np
+import numpy.typing as npt
 import pytest
 from hypothesis import given
 from hypothesis.extra.numpy import arrays
@@ -151,7 +152,7 @@ def test_invalid_value(self, dim: int) -> None:
         elements=complex_numbers(allow_infinity=False, allow_nan=False),
     ),
 )
-def test_closest_to(M: np.ndarray) -> None:
+def test_closest_to(M: npt.NDArray[np.complex128]) -> None:
     u = UnitaryMatrix.closest_to(M, (2, 2))
     assert isinstance(u, UnitaryMatrix)
 
diff --git a/tests/utils/test_math.py b/tests/utils/test_math.py
index a00a19a1f..15647aea8 100644
--- a/tests/utils/test_math.py
+++ b/tests/utils/test_math.py
@@ -4,6 +4,7 @@
 from typing import Any
 
 import numpy as np
+import numpy.typing as npt
 import pytest
 import scipy as sp
 from scipy.stats import unitary_group
@@ -17,10 +18,10 @@
 
 
 def dexpm_exact(
-    M: np.ndarray,
-    dM: np.ndarray,
+    M: npt.NDArray[np.complex128],
+    dM: npt.NDArray[np.complex128],
     term_count: int = 100,
-) -> tuple[np.ndarray, np.ndarray]:
+) -> tuple[npt.NDArray[np.complex128], npt.NDArray[np.complex128]]:
     """
     Exact matrix exponential derivative calculation.
 
@@ -45,7 +46,7 @@ class TestDexpmv:
     @pytest.mark.parametrize(
         'alpha', [np.random.random(16) for i in range(100)],
     )
-    def test_single(self, alpha: np.ndarray) -> None:
+    def test_single(self, alpha: npt.NDArray[np.float64]) -> None:
         paulis = PauliMatrices(2)
         H = paulis.dot_product(alpha)
 
@@ -59,7 +60,7 @@ def test_single(self, alpha: np.ndarray) -> None:
     @pytest.mark.parametrize(
         'alpha', [np.random.random(16) for i in range(100)],
     )
-    def test_vector(self, alpha: np.ndarray) -> None:
+    def test_vector(self, alpha: npt.NDArray[np.float64]) -> None:
         paulis = PauliMatrices(2)
         H = paulis.dot_product(alpha)
 
@@ -105,7 +106,7 @@ def test_invalid(self) -> None:
 class TestSoftmax:
 
     @pytest.mark.parametrize('x', [np.random.random(100) for i in range(100)])
-    def test_1(self, x: np.ndarray) -> None:
+    def test_1(self, x: npt.NDArray[np.float64]) -> None:
         assert np.abs(np.sum(softmax(10 * x)) - 1) < 1e-15
 
     def test_2(self) -> None:
@@ -155,7 +156,7 @@ class TestUnitaryLog:
         + PauliMatrices(4).paulis
         + [unitary_group.rvs(16) for _ in range(100)],
     )
-    def test_valid(self, reU: np.ndarray) -> None:
+    def test_valid(self, reU: npt.NDArray[np.complex128]) -> None:
         H = unitary_log_no_i(reU)
         assert np.allclose(H, H.conj().T, rtol=0, atol=1e-15)
         U = sp.linalg.expm(1j * H)
@@ -183,7 +184,7 @@ class TestPauliExpansion:
         + PauliMatrices(3).paulis
         + PauliMatrices(4).paulis,
     )
-    def test_valid(self, reH: np.ndarray) -> None:
+    def test_valid(self, reH: npt.NDArray[np.complex128]) -> None:
         alpha = pauli_expansion(reH)
         print(alpha)
         H = PauliMatrices(int(np.log2(reH.shape[0]))).dot_product(alpha)