diff --git a/examples/time_evolve_and_measure.ipynb b/examples/time_evolve_and_measure.ipynb new file mode 100644 index 000000000..acc48625d --- /dev/null +++ b/examples/time_evolve_and_measure.ipynb @@ -0,0 +1,328 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "694e0a42", + "metadata": {}, + "source": [ + "# Measure time-evolved expectation values in `qdk-chemistry`\n", + "\n", + "## This notebook demonstrates the `EvolveAndMeasure` algorithm for Hamiltonian time evolution and observable measurement using `qdk-chemistry`. It shows how to:\n", + "\n", + "1. Define a qubit Hamiltonian with time-dependent parameters\n", + "2. Build a Trotter evolution circuit\n", + "3. Run the simulation **without noise** using the QDK full-state simulator\n", + "4. Run the simulation **with noise** using both the QDK and Qiskit Aer backends\n", + "5. Optionally transpile to a target basis gate set before execution" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d7052dfd", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "from qdk_chemistry.algorithms import create\n", + "from qdk_chemistry.algorithms.time_evolution.builder.trotter import Trotter\n", + "from qdk_chemistry.algorithms.time_evolution.circuit_mapper import PauliSequenceMapper\n", + "from qdk_chemistry.data import LatticeGraph, QuantumErrorProfile, QubitHamiltonian\n", + "from qdk_chemistry.utils.model_hamiltonians import create_ising_hamiltonian\n", + "\n", + "# Reduce logging output for demo\n", + "from qdk_chemistry.utils import Logger\n", + "Logger.set_global_level(Logger.LogLevel.off)" + ] + }, + { + "cell_type": "markdown", + "id": "974e8334", + "metadata": {}, + "source": [ + "We define two qubit Hamiltonians representing alternating time-evolution steps (e.g., a driven or Floquet-like protocol). The observable is `ZZ`, measured after the full evolution sequence." + ] + }, + { + "cell_type": "markdown", + "id": "30741065", + "metadata": {}, + "source": [ + "The lists `hamiltonians` and `time_steps` define the discretized time-dependent Hamiltonian $H(t)$ as:\n", + "\n", + "$H(t_i) = H_i$, where\n", + "\n", + "$H_i \\in $ `hamiltonians` = $\\left[H_1,\\dots,H_n\\right]$,\n", + "\n", + "$t_i \\in $`time_steps` = $\\left[t_1,\\dots,t_n\\right]$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c9bdd8b5", + "metadata": {}, + "outputs": [], + "source": [ + "lattice = LatticeGraph.chain(2)\n", + "\n", + "hamiltonian_p = create_ising_hamiltonian(lattice, j=1.0, h=0.5)\n", + "hamiltonian_m = create_ising_hamiltonian(lattice, j=1.0, h=-0.5)\n", + "\n", + "steps = 20\n", + "hamiltonians = [hamiltonian_p, hamiltonian_m] * (steps // 2)\n", + "time_steps = [float((t + 1) / 10) for t in range(steps)]\n", + "\n", + "observable = QubitHamiltonian([\"ZZ\"], np.array([1.0]))" + ] + }, + { + "cell_type": "markdown", + "id": "d1affab1", + "metadata": {}, + "source": [ + "## Exact evolution via matrix exponential\n", + "\n", + "Convert the Hamiltonians to sparse matrices and compute the exact time evolution $|\\psi(t)\\rangle = \\prod_i e^{-i H_i \\Delta t_i} |\\psi_0\\rangle$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a58bb629", + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.sparse.linalg import expm_multiply\n", + "\n", + "# Convert Hamiltonians to sparse matrices\n", + "H_p = hamiltonian_p.to_matrix(sparse=True)\n", + "H_m = hamiltonian_m.to_matrix(sparse=True)\n", + "\n", + "num_qubits = H_p.shape[0].bit_length() - 1\n", + "\n", + "# Initial state |00...0⟩\n", + "psi = np.zeros(2**num_qubits, dtype=complex)\n", + "psi[0] = 1.0\n", + "\n", + "# Time-evolve: apply exp(-i H_i dt_i) for each step\n", + "dt_list = [time_steps[0]] + [time_steps[i] - time_steps[i - 1] for i in range(1, len(time_steps))]\n", + "\n", + "for ham, dt in zip(hamiltonians, dt_list):\n", + " H_sparse = H_p if ham is hamiltonian_p else H_m\n", + " psi = expm_multiply(-1j * H_sparse * dt, psi)\n", + "\n", + "# Compute exact ⟨ZZ⟩\n", + "Z = np.array([1, -1], dtype=complex)\n", + "ZZ_diag = np.kron(Z, Z)\n", + "zz_exact = np.real(np.conj(psi) @ (ZZ_diag * psi))\n", + "\n", + "print(f\"Exact ⟨ZZ⟩: {zz_exact:.6f}\")" + ] + }, + { + "cell_type": "markdown", + "id": "57064b4c", + "metadata": {}, + "source": [ + "## Noiseless simulation (QDK full-state simulator)\n", + "Set up the Trotter builder, circuit mapper, energy estimator, and the `measure_simulation` algorithm." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3ef0a94f", + "metadata": {}, + "outputs": [], + "source": [ + "evolution_builder = Trotter(num_divisions=2, order=1, optimize_term_ordering=True)\n", + "mapper = PauliSequenceMapper()\n", + "energy_estimator = create(\"energy_estimator\", \"qdk\")\n", + "measure_simulation = create(\"measure_simulation\", \"classical_sampling\")" + ] + }, + { + "cell_type": "markdown", + "id": "2ca695ec", + "metadata": {}, + "source": [ + "Run the evolution and measure `⟨ZZ⟩` without any noise using the QDK full-state simulator." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3a876d78", + "metadata": {}, + "outputs": [], + "source": [ + "circuit_executor_qdk = create(\"circuit_executor\", \"qdk_full_state_simulator\")\n", + "\n", + "measurements_noiseless = measure_simulation.run(\n", + " hamiltonians,\n", + " times=time_steps,\n", + " observables=[observable],\n", + " evolution_builder=evolution_builder,\n", + " circuit_mapper=mapper,\n", + " circuit_executor=circuit_executor_qdk,\n", + " energy_estimator=energy_estimator,\n", + " shots=10000,\n", + ")\n", + "\n", + "zz_noiseless = measurements_noiseless[0][0].energy_expectation_value\n", + "print(f\"Noiseless ⟨ZZ⟩: {zz_noiseless:.6f}\")" + ] + }, + { + "cell_type": "markdown", + "id": "6e70e127", + "metadata": {}, + "source": [ + "## Define a noise profile\n", + "\n", + "`QuantumErrorProfile` provides a backend-agnostic way to specify depolarizing noise on individual gates. This profile can be used with both the QDK and Qiskit Aer simulators.\n", + "\n", + "> **Note:** The QDK full-state simulator applies noise at the *native gate level* of the compiled QIR. If the circuit uses high-level Pauli exponentials (e.g. via `PauliSequenceMapper`), those are executed natively without decomposition into `cx`/`rz` gates — so noise on `cx` won't apply. To see noise with the QDK simulator, either use `basis_gates` to force decomposition, or define noise on the native gates (`rzz`, `rx`, etc.)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "56c5fc06", + "metadata": {}, + "outputs": [], + "source": [ + "qdk_error_profile = QuantumErrorProfile(\n", + " name=\"demo_noise\",\n", + " description=\"Light depolarizing noise on common gates\",\n", + " errors={\n", + " \"cx\": {\"type\": \"depolarizing_error\", \"rate\": 0.01, \"num_qubits\": 2},\n", + " \"cz\": {\"type\": \"depolarizing_error\", \"rate\": 0.01, \"num_qubits\": 2},\n", + " \"rzz\": {\"type\": \"depolarizing_error\", \"rate\": 0.01, \"num_qubits\": 2},\n", + " \"h\": {\"type\": \"depolarizing_error\", \"rate\": 0.001, \"num_qubits\": 1},\n", + " \"rz\": {\"type\": \"depolarizing_error\", \"rate\": 0.001, \"num_qubits\": 1},\n", + " \"rx\": {\"type\": \"depolarizing_error\", \"rate\": 0.001, \"num_qubits\": 1},\n", + " \"s\": {\"type\": \"depolarizing_error\", \"rate\": 0.001, \"num_qubits\": 1},\n", + " \"sdg\": {\"type\": \"depolarizing_error\", \"rate\": 0.001, \"num_qubits\": 1},\n", + " },\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "c27af18e", + "metadata": {}, + "source": [ + "## Noisy simulation — QDK full-state simulator\n", + "\n", + "Because `PauliSequenceMapper` produces native Pauli-exponential operations, and our noise profile includes `rzz` (the native two-qubit gate), the QDK simulator **will** apply noise to those operations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0742066a", + "metadata": {}, + "outputs": [], + "source": [ + "measurements_noisy_qdk = measure_simulation.run(\n", + " hamiltonians,\n", + " times=time_steps,\n", + " observables=[observable],\n", + " evolution_builder=evolution_builder,\n", + " circuit_mapper=mapper,\n", + " circuit_executor=circuit_executor_qdk,\n", + " energy_estimator=energy_estimator,\n", + " shots=10000,\n", + " noise=qdk_error_profile,\n", + ")\n", + "\n", + "zz_noisy_qdk = measurements_noisy_qdk[0][0].energy_expectation_value\n", + "print(f\"Noisy QDK ⟨ZZ⟩: {zz_noisy_qdk:.6f}\")" + ] + }, + { + "cell_type": "markdown", + "id": "e0518676", + "metadata": {}, + "source": [ + "## Noisy simulation — Qiskit Aer simulator\n", + "\n", + "The Qiskit Aer backend transpiles the circuit to primitive gates (`cx`, `rz`, `h`, …) before execution, so noise on those gates fires naturally." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cafb9c98", + "metadata": {}, + "outputs": [], + "source": [ + "circuit_executor_aer = create(\"circuit_executor\", \"qiskit_aer_simulator\")\n", + "\n", + "measurements_noisy_aer = measure_simulation.run(\n", + " hamiltonians,\n", + " times=time_steps,\n", + " observables=[observable],\n", + " evolution_builder=evolution_builder,\n", + " circuit_mapper=mapper,\n", + " circuit_executor=circuit_executor_aer,\n", + " energy_estimator=energy_estimator,\n", + " shots=10000,\n", + " noise=qdk_error_profile,\n", + ")\n", + "\n", + "zz_noisy_aer = measurements_noisy_aer[0][0].energy_expectation_value\n", + "print(f\"Noisy Aer ⟨ZZ⟩: {zz_noisy_aer:.6f}\")" + ] + }, + { + "cell_type": "markdown", + "id": "fdc8f4be", + "metadata": {}, + "source": [ + "## Compare results\n", + "\n", + "Side-by-side comparison of the noiseless and noisy expectation values." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8d7b777c", + "metadata": {}, + "outputs": [], + "source": [ + "print(\"Simulator ⟨ZZ⟩\")\n", + "print(\"─\" * 36)\n", + "print(f\"Exact {zz_exact: .6f}\")\n", + "print(f\"QDK (noiseless) {zz_noiseless: .6f}\")\n", + "print(f\"QDK (noisy) {zz_noisy_qdk: .6f}\")\n", + "print(f\"Aer (noisy) {zz_noisy_aer: .6f}\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "qdk_chemistry_venv", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/python/src/qdk_chemistry/algorithms/registry.py b/python/src/qdk_chemistry/algorithms/registry.py index 84e072711..079cfae20 100644 --- a/python/src/qdk_chemistry/algorithms/registry.py +++ b/python/src/qdk_chemistry/algorithms/registry.py @@ -511,11 +511,17 @@ def _register_python_factories(): from qdk_chemistry.algorithms.qubit_mapper import QubitMapperFactory # noqa: PLC0415 from qdk_chemistry.algorithms.state_preparation import StatePreparationFactory # noqa: PLC0415 from qdk_chemistry.algorithms.time_evolution.builder import TimeEvolutionBuilderFactory # noqa: PLC0415 + from qdk_chemistry.algorithms.time_evolution.circuit_mapper import ( # noqa: PLC0415 + EvolutionCircuitMapperFactory, + ) from qdk_chemistry.algorithms.time_evolution.controlled_circuit_mapper import ( # noqa: PLC0415 ControlledEvolutionCircuitMapperFactory, ) + from qdk_chemistry.algorithms.time_evolution.measure_simulation import MeasureSimulationFactory # noqa: PLC0415 register_factory(EnergyEstimatorFactory()) + register_factory(EvolutionCircuitMapperFactory()) + register_factory(MeasureSimulationFactory()) register_factory(StatePreparationFactory()) register_factory(QubitMapperFactory()) register_factory(QubitHamiltonianSolverFactory()) @@ -594,9 +600,13 @@ def _register_python_algorithms(): from qdk_chemistry.algorithms.time_evolution.builder.trotter import ( # noqa: PLC0415 Trotter, ) + from qdk_chemistry.algorithms.time_evolution.circuit_mapper import ( # noqa: PLC0415 + PauliSequenceMapper as EvolutionPauliSequenceMapper, + ) from qdk_chemistry.algorithms.time_evolution.controlled_circuit_mapper import ( # noqa: PLC0415 - PauliSequenceMapper, + PauliSequenceMapper as ControlledPauliSequenceMapper, ) + from qdk_chemistry.algorithms.time_evolution.measure_simulation import EvolveAndMeasure # noqa: PLC0415 register(lambda: QdkEnergyEstimator()) register(lambda: SparseIsometryGF2XStatePreparation()) @@ -606,7 +616,9 @@ def _register_python_algorithms(): register(lambda: Trotter()) register(lambda: QDrift()) register(lambda: PartiallyRandomized()) - register(lambda: PauliSequenceMapper()) + register(lambda: EvolutionPauliSequenceMapper()) + register(lambda: ControlledPauliSequenceMapper()) + register(lambda: EvolveAndMeasure()) register(lambda: QdkFullStateSimulator()) register(lambda: QdkSparseStateSimulator()) register(lambda: IterativePhaseEstimation()) diff --git a/python/src/qdk_chemistry/algorithms/time_evolution/builder/trotter.py b/python/src/qdk_chemistry/algorithms/time_evolution/builder/trotter.py index 40839bb26..f1e1dcd5d 100644 --- a/python/src/qdk_chemistry/algorithms/time_evolution/builder/trotter.py +++ b/python/src/qdk_chemistry/algorithms/time_evolution/builder/trotter.py @@ -18,6 +18,8 @@ # Licensed under the MIT License. See LICENSE.txt in the project root for license information. # -------------------------------------------------------------------------------------------- +import numpy as np + from qdk_chemistry.algorithms.time_evolution.builder.base import TimeEvolutionBuilder from qdk_chemistry.algorithms.time_evolution.builder.trotter_error import ( trotter_steps_commutator, @@ -45,6 +47,7 @@ def __init__(self): num_divisions: Explicit number of divisions within a Trotter step (0 means automatic). error_bound: Strategy for computing the Trotter error bound ("commutator" or "naive"). weight_threshold: The absolute threshold for filtering small coefficients. + optimize_term_ordering: Whether to group commuting terms and execute them in parallel. """ super().__init__() @@ -71,6 +74,12 @@ def __init__(self): self._set_default( "weight_threshold", "float", 1e-12, "The absolute threshold for filtering small coefficients." ) + self._set_default( + "optimize_term_ordering", + "bool", + False, + "Whether to group commuting terms and execute them in parallel.", + ) class Trotter(TimeEvolutionBuilder): @@ -84,6 +93,7 @@ def __init__( num_divisions: int = 0, error_bound: str = "commutator", weight_threshold: float = 1e-12, + optimize_term_ordering: bool = False, ): r"""Initialize Trotter builder with specified Trotter decomposition settings. @@ -131,6 +141,7 @@ def __init__( (default, tighter) or ``"naive"``. weight_threshold: Absolute threshold for filtering small Hamiltonian coefficients. Defaults to 1e-12. + optimize_term_ordering: Whether to group commuting terms and execute them in parallel. Defaults to False. """ super().__init__() @@ -140,6 +151,7 @@ def __init__( self._settings.set("num_divisions", num_divisions) self._settings.set("error_bound", error_bound) self._settings.set("weight_threshold", weight_threshold) + self._settings.set("optimize_term_ordering", optimize_term_ordering) def _run_impl(self, qubit_hamiltonian: QubitHamiltonian, time: float) -> TimeEvolutionUnitary: """Construct the time evolution unitary using Trotter decomposition. @@ -187,7 +199,11 @@ def _trotter(self, qubit_hamiltonian: QubitHamiltonian, time: float) -> TimeEvol delta = time / num_divisions - terms = self._decompose_trotter_step(qubit_hamiltonian, time=delta, atol=weight_threshold) + optimize_term_ordering = self._settings.get("optimize_term_ordering") + + terms = self._decompose_trotter_step( + qubit_hamiltonian, time=delta, atol=weight_threshold, optimize_term_ordering=optimize_term_ordering + ) num_qubits = qubit_hamiltonian.num_qubits @@ -237,7 +253,12 @@ def _resolve_num_divisions(self, qubit_hamiltonian: QubitHamiltonian, time: floa return max(manual, auto) def _decompose_trotter_step( - self, qubit_hamiltonian: QubitHamiltonian, time: float, *, atol: float = 1e-12 + self, + qubit_hamiltonian: QubitHamiltonian, + time: float, + *, + atol: float = 1e-12, + optimize_term_ordering: bool = False, ) -> list[ExponentiatedPauliTerm]: """Decompose a single Trotter step into exponentiated Pauli terms. @@ -249,6 +270,8 @@ def _decompose_trotter_step( time: The evolution time for the single step. atol: Absolute tolerance for filtering small coefficients. + optimize_term_ordering: Whether to group commuting terms together + and further subgroup into parallelizable layers. Returns: A list of ``ExponentiatedPauliTerm`` representing the decomposed terms. @@ -259,8 +282,6 @@ def _decompose_trotter_step( if not qubit_hamiltonian.is_hermitian(tolerance=atol): raise ValueError("Non-Hermitian Hamiltonian: coefficients have nonzero imaginary parts.") - order = self._settings.get("order") - coeffs = list(qubit_hamiltonian.get_real_coefficients(tolerance=atol)) # If there are no coefficients (e.g., empty Hamiltonian or all filtered by atol), # there is nothing to decompose; return the empty list of terms. @@ -268,32 +289,49 @@ def _decompose_trotter_step( Logger.warn("No coefficients above the tolerance; returning empty term list.") return terms + order = self._settings.get("order") + grouped_hamiltonians = self._group_terms(qubit_hamiltonian, optimize_term_ordering=optimize_term_ordering) + if order == 1: - for label, coeff in coeffs: - mapping = self._pauli_label_to_map(label) - angle = coeff * time - terms.append(ExponentiatedPauliTerm(pauli_term=mapping, angle=angle)) + for group in grouped_hamiltonians: + for subgroup in group: + terms.extend( + self._exponentiate_commuting( + subgroup, + time=time, + atol=atol, + ) + ) + # order = 2 or order = 2k with k>1 else: - # \prod_{i=1}^{L-1} e^{-iH_i t/(2n)} - for label, coeff in coeffs[:-1]: - mapping = self._pauli_label_to_map(label) - angle = coeff * time / 2 - terms.append(ExponentiatedPauliTerm(pauli_term=mapping, angle=angle)) - # e^{-iH_L t/n} - label, coeff = coeffs[-1] - mapping = self._pauli_label_to_map(label) - angle = coeff * time - terms.append(ExponentiatedPauliTerm(pauli_term=mapping, angle=angle)) + # \prod e^{-iH_i t/(2n)} for all Hamiltonian terms groups except the last one + # which is executed in the middle as e^{-iH_i t/n} + terms_without_last_group = [] + for group in grouped_hamiltonians[:-1]: + for subgroup in group: + terms_without_last_group.extend( + self._exponentiate_commuting( + subgroup, + time=time / 2, + atol=atol, + ) + ) - # \prod_{i=L-1}^1 e^{-iH_i t/(2n)} - for label, coeff in reversed(coeffs[:-1]): - mapping = self._pauli_label_to_map(label) - angle = coeff * time / 2 - terms.append(ExponentiatedPauliTerm(pauli_term=mapping, angle=angle)) + terms.extend(terms_without_last_group) + # e^{-iH_i t/n} for all terms in the last group + for subgroup in grouped_hamiltonians[-1]: + terms.extend( + self._exponentiate_commuting( + subgroup, + time=time, + atol=atol, + ) + ) + terms.extend(terms_without_last_group[::-1]) # reverse all but the last group and append # Construct order 2k formula bottom up dynamic-programming style - if order > 2: + if order > 2 and order % 2 == 0: step_terms = terms.copy() for k in range(2, int(order / 2) + 1): u_k = 1 / (4 - 4 ** (1 / (2 * k - 1))) @@ -342,6 +380,137 @@ def _decompose_trotter_step( else: merged_terms.append(term) terms = merged_terms + + return terms + + def _group_terms( + self, + qubit_hamiltonian: QubitHamiltonian, + *, + optimize_term_ordering: bool = True, + ) -> list[list[QubitHamiltonian]]: + """Group Hamiltonian terms into commuting and concurrent sets. + + When *optimize_term_ordering* is ``True``: + 1. Partition using :meth:`QubitHamiltonian.group_commuting`. + 2. Within each group, merge terms with identical Pauli strings + by summing their coefficients. + 3. Within each group, split terms into parallelizable layers + (terms whose Pauli supports are disjoint). + 4. Move the group with the most multi-qubit terms to the end. + + When ``False``, each Pauli string is returned as its own + single-term group with no reordering. + + Args: + qubit_hamiltonian: The qubit Hamiltonian to group. + optimize_term_ordering: Whether to group commuting terms together. + + Returns: + A list of groups, where each group is a list of + ``QubitHamiltonian`` sub-groups (parallelizable layers). + + """ + if not optimize_term_ordering: + return [ + [ + QubitHamiltonian( + pauli_strings=[label], + coefficients=[coeff], + encoding=qubit_hamiltonian.encoding, + ) + ] + for label, coeff in zip(qubit_hamiltonian.pauli_strings, qubit_hamiltonian.coefficients, strict=True) + ] + + # Sort terms so that Pauli strings acting on more qubits appear first. + num_non_identity = [sum(c != "I" for c in ps) for ps in qubit_hamiltonian.pauli_strings] + sorted_indices = sorted(range(len(num_non_identity)), key=lambda i: num_non_identity[i], reverse=True) + qubit_hamiltonian = QubitHamiltonian( + pauli_strings=[qubit_hamiltonian.pauli_strings[i] for i in sorted_indices], + coefficients=np.asarray([qubit_hamiltonian.coefficients[i] for i in sorted_indices]), + encoding=qubit_hamiltonian.encoding, + ) + + sub_hamiltonians = qubit_hamiltonian.group_commuting(qubit_wise=False) + + result: list[list[QubitHamiltonian]] = [] + for sub_h in sub_hamiltonians: + # Merge terms with identical Pauli strings. + merged: dict[str, complex] = {} + for label, coeff in zip(sub_h.pauli_strings, sub_h.coefficients, strict=True): + merged[label] = merged.get(label, 0.0) + coeff + labels = list(merged.keys()) + coeffs = list(merged.values()) + + # Split into parallelizable layers (disjoint qubit supports). + # Each layer becomes its own sub-group consisting of terms whose + # supports are mutually disjoint, allowing them to be applied in parallel. + pauli_maps = [self._pauli_label_to_map(label) for label in labels] + layers_indices: list[list[int]] = [] + layers_occupied: list[set[int]] = [] + for i, pm in enumerate(pauli_maps): + qubits = set(pm.keys()) + placed = False + for layer_i, layer_occ in enumerate(layers_occupied): + if qubits.isdisjoint(layer_occ): + layers_indices[layer_i].append(i) + layer_occ.update(qubits) + placed = True + break + if not placed: + layers_indices.append([i]) + layers_occupied.append(set(qubits)) + + outer_group: list[QubitHamiltonian] = [] + for layer in layers_indices: + outer_group.append( + QubitHamiltonian( + pauli_strings=[labels[i] for i in layer], + coefficients=np.asarray([coeffs[i] for i in layer]), + encoding=sub_h.encoding, + ) + ) + result.append(outer_group) + + # Move the group with the most multi-qubit terms to the end. + def _multi_qubit_count(group: list[QubitHamiltonian]) -> int: + return sum(1 for sub_h in group for label in sub_h.pauli_strings if sum(c != "I" for c in label) > 1) + + max_idx = max(range(len(result)), key=lambda i: _multi_qubit_count(result[i])) + result.append(result.pop(max_idx)) + + return result + + def _exponentiate_commuting( + self, + group: QubitHamiltonian, + time: float, + *, + atol: float = 1e-12, + ) -> list[ExponentiatedPauliTerm]: + r"""Exponentiate a group of commuting Pauli terms. + + Each term :math:`P_j` with coefficient :math:`c_j` is converted to + the rotation :math:`e^{-i\,c_j\,t\,P_j}`. Because all terms in the + group commute and :meth:`_group_terms` ensures they have disjoint + qubit supports, the rotations can be applied in any order. + + Args: + group: The group of commuting Hamiltonian terms to exponentiate. + time: The evolution time used to compute rotation angles + (:math:`\theta_j = c_j \cdot t`). + atol: Absolute tolerance for filtering small coefficients. + + Returns: + A flat list of :class:`ExponentiatedPauliTerm`. + + """ + terms: list[ExponentiatedPauliTerm] = [] + for label, coeff in group.get_real_coefficients(tolerance=atol): + mapping = self._pauli_label_to_map(label) + angle = coeff * time + terms.append(ExponentiatedPauliTerm(pauli_term=mapping, angle=angle)) return terms def name(self) -> str: diff --git a/python/src/qdk_chemistry/algorithms/time_evolution/circuit_mapper/__init__.py b/python/src/qdk_chemistry/algorithms/time_evolution/circuit_mapper/__init__.py new file mode 100644 index 000000000..59343771f --- /dev/null +++ b/python/src/qdk_chemistry/algorithms/time_evolution/circuit_mapper/__init__.py @@ -0,0 +1,15 @@ +"""QDK/Chemistry time evolution circuit mapper module.""" + +# -------------------------------------------------------------------------------------------- +# Copyright (c) Microsoft Corporation. All rights reserved. +# Licensed under the MIT License. See LICENSE.txt in the project root for license information. +# -------------------------------------------------------------------------------------------- + +from .base import EvolutionCircuitMapperFactory +from .pauli_sequence_mapper import PauliSequenceMapper, PauliSequenceMapperSettings + +__all__ = [ + "EvolutionCircuitMapperFactory", + "PauliSequenceMapper", + "PauliSequenceMapperSettings", +] diff --git a/python/src/qdk_chemistry/algorithms/time_evolution/circuit_mapper/base.py b/python/src/qdk_chemistry/algorithms/time_evolution/circuit_mapper/base.py new file mode 100644 index 000000000..70d299ad0 --- /dev/null +++ b/python/src/qdk_chemistry/algorithms/time_evolution/circuit_mapper/base.py @@ -0,0 +1,49 @@ +"""QDK/Chemistry time evolution unitary circuit mapper abstractions.""" + +# -------------------------------------------------------------------------------------------- +# Copyright (c) Microsoft Corporation. All rights reserved. +# Licensed under the MIT License. See LICENSE.txt in the project root for license information. +# -------------------------------------------------------------------------------------------- + +from abc import abstractmethod + +from qdk_chemistry.algorithms.base import Algorithm, AlgorithmFactory +from qdk_chemistry.data import Circuit +from qdk_chemistry.data.time_evolution.base import TimeEvolutionUnitary + +__all__: list[str] = ["EvolutionCircuitMapper", "EvolutionCircuitMapperFactory"] + + +class EvolutionCircuitMapper(Algorithm): + """Base class for time evolution circuit mapper in QDK/Chemistry algorithms.""" + + def __init__(self): + """Initialize the EvolutionCircuitMapper.""" + super().__init__() + + @abstractmethod + def _run_impl(self, evolution: TimeEvolutionUnitary, *args, **kwargs) -> Circuit: + """Construct a Circuit representing the unitary for the given TimeEvolutionUnitary. + + Args: + evolution: The time evolution unitary. + *args: Positional arguments, where the first argument is expected to be the + time evolution unitary. + **kwargs: Additional keyword arguments for concrete implementation. + + Returns: + Circuit: A Circuit representing the unitary for the given TimeEvolutionUnitary. + + """ + + +class EvolutionCircuitMapperFactory(AlgorithmFactory): + """Factory class for creating EvolutionCircuitMapper instances.""" + + def algorithm_type_name(self) -> str: + """Return evolution_circuit_mapper as the algorithm type name.""" + return "evolution_circuit_mapper" + + def default_algorithm_name(self) -> str: + """Return pauli_sequence as the default algorithm name.""" + return "pauli_sequence" diff --git a/python/src/qdk_chemistry/algorithms/time_evolution/circuit_mapper/pauli_sequence_mapper.py b/python/src/qdk_chemistry/algorithms/time_evolution/circuit_mapper/pauli_sequence_mapper.py new file mode 100644 index 000000000..f2056e0f5 --- /dev/null +++ b/python/src/qdk_chemistry/algorithms/time_evolution/circuit_mapper/pauli_sequence_mapper.py @@ -0,0 +1,115 @@ +"""QDK/Chemistry sequence structure evolution circuit mapper.""" + +# -------------------------------------------------------------------------------------------- +# Copyright (c) Microsoft Corporation. All rights reserved. +# Licensed under the MIT License. See LICENSE.txt in the project root for license information. +# -------------------------------------------------------------------------------------------- + +from qdk import qsharp + +from qdk_chemistry.data import Settings +from qdk_chemistry.data.circuit import Circuit +from qdk_chemistry.data.time_evolution.base import TimeEvolutionUnitary +from qdk_chemistry.data.time_evolution.containers.pauli_product_formula import ( + PauliProductFormulaContainer, +) +from qdk_chemistry.utils.qsharp import QSHARP_UTILS + +from .base import EvolutionCircuitMapper + +__all__: list[str] = ["PauliSequenceMapper", "PauliSequenceMapperSettings"] + + +class PauliSequenceMapperSettings(Settings): + """Settings for PauliSequenceMapper.""" + + def __init__(self): + """Initialize PauliSequenceMapperSettings with default values.""" + super().__init__() + + +class PauliSequenceMapper(EvolutionCircuitMapper): + r"""Evolution circuit mapper using Pauli product formula term sequences. + + Given a time-evolution operator expressed as a Pauli product formula + :math:`U(t) \approx \left[ U_{\mathrm{step}}(t / r) \right]^{r}`, this mapper constructs + a :math:`U(t)` using the following pattern: + + 1. Each Pauli operator :math:`P_j` is basis-rotated into the :math:`Z` basis. + 2. Qubits involved in :math:`P_j` are entangled into a sequence using CNOT gates. + 3. A :math:`R_z` rotation implements + :math:`e^{-i\,\theta_j\,P_j} \;\rightarrow\; R_z(2 \theta_j)`. + 4. The basis rotations and entangling operations are uncomputed. + + Notes: + * Requires a ``PauliProductFormulaContainer`` for the time evolution unitary. + + """ + + def __init__(self): + """Initialize the PauliSequenceMapper.""" + super().__init__() + self._settings = PauliSequenceMapperSettings() + + def name(self) -> str: + """Return the algorithm name.""" + return "pauli_sequence" + + def type_name(self) -> str: + """Return evolution_circuit_mapper as the algorithm type name.""" + return "evolution_circuit_mapper" + + def _run_impl(self, evolution: TimeEvolutionUnitary) -> Circuit: + r"""Construct a quantum circuit implementing the time evolution unitary. + + Args: + evolution: The time evolution unitary containing the Hamiltonian + and evolution parameters. + + Returns: + Circuit: A quantum circuit implementing the unitary :math:`U` + where :math:`U` is the time evolution operator :math:`\exp(-i H t)`. + + Raises: + ValueError: If the time evolution unitary container type is not supported. + + """ + unitary_container = evolution.get_container() + if not isinstance(unitary_container, PauliProductFormulaContainer): + raise ValueError( + f"The {evolution.get_container_type()} container type is not supported. " + "PauliSequenceMapper only supports PauliProductFormula container for time evolution unitary." + ) + + pauli_terms: list[list[qsharp.Pauli]] = [] + angles: list[float] = [] + for term in unitary_container.step_terms: + base_terms = [qsharp.Pauli.I] * unitary_container.num_qubits + for index, pauli in term.pauli_term.items(): + base_terms[index] = getattr(qsharp.Pauli, pauli) + pauli_terms.append(base_terms.copy()) + angles.append(term.angle) + + evo_params = { + "pauliExponents": pauli_terms, + "pauliCoefficients": angles, + "repetitions": unitary_container.step_reps, + } + + target_indices = list(range(unitary_container.num_qubits)) + + qsc = qsharp.circuit( + QSHARP_UTILS.PauliExp.MakeRepPauliExpCircuit, + evo_params, + target_indices, + ) + + qir = qsharp.compile( + QSHARP_UTILS.PauliExp.MakeRepPauliExpCircuit, + evo_params, + target_indices, + ) + + evolution_op = QSHARP_UTILS.PauliExp.MakeRepPauliExpOp(evo_params) + + return Circuit(qsharp=qsc, qir=qir, qsharp_op=evolution_op) diff --git a/python/src/qdk_chemistry/algorithms/time_evolution/measure_simulation/__init__.py b/python/src/qdk_chemistry/algorithms/time_evolution/measure_simulation/__init__.py new file mode 100644 index 000000000..eab2c9b4f --- /dev/null +++ b/python/src/qdk_chemistry/algorithms/time_evolution/measure_simulation/__init__.py @@ -0,0 +1,10 @@ +"""QDK/Chemistry evolve-and-measure algorithms module.""" + +# -------------------------------------------------------------------------------------------- +# Copyright (c) Microsoft Corporation. All rights reserved. +# Licensed under the MIT License. See LICENSE.txt in the project root for license information. +# -------------------------------------------------------------------------------------------- +from .base import MeasureSimulationFactory +from .evolve_and_measure import EvolveAndMeasure + +__all__: list[str] = ["EvolveAndMeasure", "MeasureSimulationFactory"] diff --git a/python/src/qdk_chemistry/algorithms/time_evolution/measure_simulation/base.py b/python/src/qdk_chemistry/algorithms/time_evolution/measure_simulation/base.py new file mode 100644 index 000000000..1f64cbc2c --- /dev/null +++ b/python/src/qdk_chemistry/algorithms/time_evolution/measure_simulation/base.py @@ -0,0 +1,242 @@ +"""QDK/Chemistry evolve-and-measure abstractions and utilities.""" + +# -------------------------------------------------------------------------------------------- +# Copyright (c) Microsoft Corporation. All rights reserved. +# Licensed under the MIT License. See LICENSE.txt in the project root for license information. +# -------------------------------------------------------------------------------------------- + +import re +from abc import abstractmethod + +from qdk import qsharp + +from qdk_chemistry.algorithms.base import Algorithm, AlgorithmFactory +from qdk_chemistry.algorithms.circuit_executor.base import CircuitExecutor +from qdk_chemistry.algorithms.energy_estimator.energy_estimator import EnergyEstimator +from qdk_chemistry.algorithms.time_evolution.builder.base import TimeEvolutionBuilder +from qdk_chemistry.algorithms.time_evolution.circuit_mapper.base import EvolutionCircuitMapper +from qdk_chemistry.data import ( + Circuit, + EnergyExpectationResult, + MeasurementData, + QuantumErrorProfile, + QubitHamiltonian, + Settings, + TimeEvolutionUnitary, +) +from qdk_chemistry.utils.qsharp import QSHARP_UTILS + +__all__: list[str] = ["MeasureSimulation", "MeasureSimulationFactory", "MeasureSimulationSettings"] + +_QASM_QUBIT_DECLARATION_PATTERN = re.compile(r"\bqubit\s*\[(\d+)\]") +_QIR_REQUIRED_NUM_QUBITS_PATTERN = re.compile(r'"required_num_qubits"="(\d+)"') + + +class MeasureSimulationSettings(Settings): + """Settings for evolve-and-measure simulation.""" + + def __init__(self): + """Initialize defaults for evolve-and-measure simulation.""" + super().__init__() + + +class MeasureSimulation(Algorithm): + """Abstract base class for Hamiltonian evolution and observable measurement.""" + + def __init__(self): + """Initialize the evolve-and-measure simulation settings.""" + super().__init__() + self._settings = MeasureSimulationSettings() + + def type_name(self) -> str: + """Return the algorithm type name as measure_simulation.""" + return "measure_simulation" + + @abstractmethod + def _run_impl( + self, + qubit_hamiltonians: list[QubitHamiltonian], + times: list[float], + observables: list[QubitHamiltonian], + *, + state_prep: Circuit | None = None, + evolution_builder: TimeEvolutionBuilder, + circuit_mapper: EvolutionCircuitMapper, + shots: int = 1000, + circuit_executor: CircuitExecutor, + energy_estimator: EnergyEstimator, + noise: QuantumErrorProfile | None = None, + basis_gates: list[str] | None = None, + ) -> list[tuple[EnergyExpectationResult, MeasurementData]]: + """Run evolve-and-measure simulation. + + Args: + qubit_hamiltonians: List of Hamiltonians used to build time evolution. + times: List of times to evolve under the Hamiltonians. + observables: List of observable Hamiltonians to measure after evolution. + state_prep: Optional circuit that prepares the initial state before time evolution. + evolution_builder: Time-evolution builder. + circuit_mapper: Mapper for time-evolution unitary to circuit. + shots: Number of shots to use for measurement. + circuit_executor: Circuit executor backend. + energy_estimator: Energy estimator algorithm. + noise: Optional noise profile. + basis_gates: Optional list of basis gates to transpile the circuit into before execution. + + Returns: + A list of tuples containing ``EnergyExpectationResult`` and ``MeasurementData`` objects. + + """ + + def _create_time_evolution( + self, + qubit_hamiltonian: QubitHamiltonian, + time: float, + evolution_builder: TimeEvolutionBuilder, + ) -> TimeEvolutionUnitary: + """Create the time-evolution unitary for current settings.""" + return evolution_builder.run(qubit_hamiltonian, time) + + def _map_time_evolution_to_circuit( + self, + evolution: TimeEvolutionUnitary, + circuit_mapper: EvolutionCircuitMapper, + ) -> Circuit: + """Map a time-evolution unitary into an executable circuit.""" + return circuit_mapper.run(evolution) + + def _prepend_state_prep_circuit(self, state_prep: Circuit, circuit: Circuit) -> Circuit: + state_prep_op = getattr(state_prep, "_qsharp_op", None) + circuit_op = getattr(circuit, "_qsharp_op", None) + if state_prep_op is None or circuit_op is None: + raise RuntimeError("State-preparation circuit composition requires Q# operations on both circuits.") + + if state_prep.encoding is not None and circuit.encoding is not None and state_prep.encoding != circuit.encoding: + raise ValueError( + "State-preparation circuit and evolution circuit use different encodings " + f"('{state_prep.encoding}' and '{circuit.encoding}')." + ) + + num_qubits = None + for representation in (state_prep.qir, circuit.qir, state_prep.qasm, circuit.qasm): + if representation is None: + continue + + match = _QIR_REQUIRED_NUM_QUBITS_PATTERN.search(str(representation)) + if match: + candidate = int(match.group(1)) + else: + qasm_match = _QASM_QUBIT_DECLARATION_PATTERN.search(str(representation)) + if not qasm_match: + continue + candidate = int(qasm_match.group(1)) + + if num_qubits is None: + num_qubits = candidate + elif num_qubits != candidate: + raise ValueError( + "State-preparation circuit and evolution circuit must act on the same number of qubits " + f"(received {num_qubits} and {candidate})." + ) + + if num_qubits is None: + raise RuntimeError("Unable to infer the number of qubits needed to compose the Q# circuits.") + + target_indices = list(range(num_qubits)) + combined_encoding = circuit.encoding if circuit.encoding is not None else state_prep.encoding + qsharp_circuit = qsharp.circuit( + QSHARP_UTILS.CircuitComposition.MakeSequentialCircuit, + state_prep_op, + circuit_op, + target_indices, + ) + qir = qsharp.compile( + QSHARP_UTILS.CircuitComposition.MakeSequentialCircuit, + state_prep_op, + circuit_op, + target_indices, + ) + + return Circuit( + qsharp=qsharp_circuit, + qir=qir, + qsharp_op=QSHARP_UTILS.CircuitComposition.MakeSequentialOp(state_prep_op, circuit_op), + encoding=combined_encoding, + ) + + @staticmethod + def _transpile_to_basis_gates(circuit: Circuit, basis_gates: list[str]) -> Circuit: + """Transpile a Circuit to a target basis gate set using the qdk-chemistry transpiler. + + Args: + circuit: The circuit to transpile. + basis_gates: Target basis gates (e.g. ``["cx", "rz", "h", "x"]``). + + Returns: + A new ``Circuit`` restricted to the requested basis gates. + + """ + try: + from qiskit import qasm3, transpile # noqa: PLC0415 + from qiskit.transpiler import PassManager # noqa: PLC0415 + + from qdk_chemistry.plugins.qiskit._interop.transpiler import ( # noqa: PLC0415 + MergeZBasisRotations, + RemoveZBasisOnZeroState, + SubstituteCliffordRz, + ) + except ImportError as exc: + raise RuntimeError( + "Qiskit is required to transpile circuits to the requested basis_gates, " + "but it is not installed. Please install the 'qiskit' package to use " + "the basis_gates option." + ) from exc + + qc = circuit.get_qiskit_circuit() + qc = transpile(qc, basis_gates=basis_gates, optimization_level=3) + + pm = PassManager( + [ + MergeZBasisRotations(), + SubstituteCliffordRz(), + RemoveZBasisOnZeroState(), + ] + ) + qc = pm.run(qc) + + return Circuit(qasm=qasm3.dumps(qc), encoding=circuit.encoding) + + def _measure_observable( + self, + circuit: Circuit, + observable: QubitHamiltonian, + circuit_executor: CircuitExecutor, + energy_estimator: EnergyEstimator, + shots: int = 1000, + noise: QuantumErrorProfile | None = None, + ) -> tuple[EnergyExpectationResult, MeasurementData]: + """Measure a qubit observable on the provided circuit state.""" + energy_result, measurement_data = energy_estimator.run( + circuit, + observable, + circuit_executor, + total_shots=shots, + noise_model=noise, + ) + return energy_result, measurement_data + + +class MeasureSimulationFactory(AlgorithmFactory): + """Factory class for creating evolve-and-measure algorithm instances.""" + + def __init__(self): + """Initialize the MeasureSimulationFactory.""" + super().__init__() + + def algorithm_type_name(self) -> str: + """Return the algorithm type name as measure_simulation.""" + return "measure_simulation" + + def default_algorithm_name(self) -> str: + """Return classical sampling as the default algorithm name.""" + return "classical_sampling" diff --git a/python/src/qdk_chemistry/algorithms/time_evolution/measure_simulation/evolve_and_measure.py b/python/src/qdk_chemistry/algorithms/time_evolution/measure_simulation/evolve_and_measure.py new file mode 100644 index 000000000..ae48b9ad3 --- /dev/null +++ b/python/src/qdk_chemistry/algorithms/time_evolution/measure_simulation/evolve_and_measure.py @@ -0,0 +1,136 @@ +"""Hamiltonian evolution + observable measurement implementation.""" + +# -------------------------------------------------------------------------------------------- +# Copyright (c) Microsoft Corporation. All rights reserved. +# Licensed under the MIT License. See LICENSE.txt in the project root for license information. +# -------------------------------------------------------------------------------------------- + +from qdk_chemistry.algorithms.circuit_executor.base import CircuitExecutor +from qdk_chemistry.algorithms.energy_estimator.energy_estimator import EnergyEstimator +from qdk_chemistry.algorithms.time_evolution.builder.base import TimeEvolutionBuilder +from qdk_chemistry.algorithms.time_evolution.circuit_mapper.base import EvolutionCircuitMapper +from qdk_chemistry.data import ( + Circuit, + EnergyExpectationResult, + MeasurementData, + QuantumErrorProfile, + QubitHamiltonian, + TimeEvolutionUnitary, +) +from qdk_chemistry.utils import Logger + +from .base import MeasureSimulation, MeasureSimulationSettings + +__all__: list[str] = ["EvolveAndMeasure", "EvolveAndMeasureSettings"] + + +class EvolveAndMeasureSettings(MeasureSimulationSettings): + """Settings for the EvolveAndMeasure algorithm.""" + + def __init__(self): + """Initialize the settings for EvolveAndMeasure.""" + super().__init__() + + +class EvolveAndMeasure(MeasureSimulation): + """Evolve under a Hamiltonian and measure a target observable.""" + + def __init__(self): + """Initialize EvolveAndMeasure with the given settings.""" + Logger.trace_entering() + super().__init__() + self._settings = EvolveAndMeasureSettings() + + def _run_impl( + self, + qubit_hamiltonians: list[QubitHamiltonian], + times: list[float], + observables: list[QubitHamiltonian], + *, + state_prep: Circuit | None = None, + evolution_builder: TimeEvolutionBuilder, + circuit_mapper: EvolutionCircuitMapper, + shots: int = 1000, + circuit_executor: CircuitExecutor, + energy_estimator: EnergyEstimator, + noise: QuantumErrorProfile | None = None, + basis_gates: list[str] | None = None, + ) -> list[tuple[EnergyExpectationResult, MeasurementData]]: + """Run evolve-and-measure simulation. + + Args: + qubit_hamiltonians: List of Hamiltonians used to build time evolution. + times: Monotonically-increasing list of times to evolve under the Hamiltonians. + observables: List of observable Hamiltonians to measure after evolution. + state_prep: Optional circuit that prepares the initial state before time evolution. + evolution_builder: Time-evolution builder. + circuit_mapper: Mapper for time-evolution unitary to circuit. + shots: Number of shots to use for measurement. + circuit_executor: Circuit executor backend. + energy_estimator: Energy estimator algorithm. + noise: Optional noise profile. + basis_gates: Optional list of basis gates to transpile the circuit into before execution. + + Returns: + A list of tuples containing ``EnergyExpectationResult`` and ``MeasurementData`` objects. + + Raises: + ValueError: If ``qubit_hamiltonians`` is empty. + + """ + if not qubit_hamiltonians: + raise ValueError("qubit_hamiltonians must not be empty.") + if not times: + raise ValueError("times must not be empty.") + if len(qubit_hamiltonians) != len(times): + raise ValueError("qubit_hamiltonians and times must have the same length.") + if times != sorted(times): + raise ValueError("times must be monotonically increasing.") + + # Ensure all Hamiltonians and observables have the same number of qubits. + reference_num_qubits = qubit_hamiltonians[0].num_qubits + for hamiltonian in qubit_hamiltonians[1:]: + if hamiltonian.num_qubits != reference_num_qubits: + raise ValueError("All Hamiltonians and observables must have the same number of qubits.") + for observable in observables: + if observable.num_qubits != reference_num_qubits: + raise ValueError("All Hamiltonians and observables must have the same number of qubits.") + evolution = self._create_time_evolution(qubit_hamiltonians[0], times[0], evolution_builder) + + for i in range(1, len(qubit_hamiltonians)): + qubit_hamiltonian = qubit_hamiltonians[i] + time = times[i] + delta_t = time - times[i - 1] + + evolution = TimeEvolutionUnitary( + evolution.get_container().combine( + self._create_time_evolution(qubit_hamiltonian, delta_t, evolution_builder).get_container(), + ) + ) + + evolution_circuit = self._map_time_evolution_to_circuit(evolution, circuit_mapper) + + if state_prep is not None: + evolution_circuit = self._prepend_state_prep_circuit(state_prep, evolution_circuit) + + # Transpile to basis gates + if basis_gates is not None: + evolution_circuit = self._transpile_to_basis_gates(evolution_circuit, basis_gates) + + measurements = [] + for observable in observables: + measurements.append( + self._measure_observable( + circuit=evolution_circuit, + shots=shots, + observable=observable, + circuit_executor=circuit_executor, + energy_estimator=energy_estimator, + noise=noise, + ) + ) + return measurements + + def name(self) -> str: + """Return ``classical_sampling`` as the algorithm name.""" + return "classical_sampling" diff --git a/python/src/qdk_chemistry/data/time_evolution/containers/pauli_product_formula.py b/python/src/qdk_chemistry/data/time_evolution/containers/pauli_product_formula.py index 3d7de3eb5..f6850ef1f 100644 --- a/python/src/qdk_chemistry/data/time_evolution/containers/pauli_product_formula.py +++ b/python/src/qdk_chemistry/data/time_evolution/containers/pauli_product_formula.py @@ -125,6 +125,51 @@ def reorder_terms(self, permutation: list[int]) -> "PauliProductFormulaContainer num_qubits=self._num_qubits, ) + def combine(self, other_container: "PauliProductFormulaContainer", atol=1e-12) -> "PauliProductFormulaContainer": + """Combine two Trotter evolutions, merging adjacent identical Pauli terms. + + The terms from ``self`` (repeated ``step_reps`` times) are followed by the + terms from ``other_container`` (also repeated according to its + ``step_reps``). When two consecutive terms act with the same Pauli operator + string (i.e., have identical ``pauli_term`` dictionaries), their rotation + angles are summed into a single ``ExponentiatedPauliTerm``. If the summed + angle has magnitude less than ``atol``, the resulting term is removed. + + Args: + other_container: The second ``PauliProductFormulaContainer`` appended + after this container. + atol: Absolute tolerance used when deciding whether a merged term with + a small rotation angle should be dropped. + + Returns: + A single ``PauliProductFormulaContainer`` representing the combined + evolution with adjacent identical terms fused. + + """ + if self.num_qubits != other_container.num_qubits: + raise ValueError( + f"Cannot combine PauliProductFormulaContainer instances with different " + f"num_qubits (self.num_qubits={self.num_qubits}, " + f"other_container.num_qubits={other_container.num_qubits})." + ) + + merged: list[ExponentiatedPauliTerm] = [] + for step_terms, step_reps in ( + (self.step_terms, self.step_reps), + (other_container.step_terms, other_container.step_reps), + ): + for _ in range(step_reps): + for term in step_terms: + if merged and merged[-1].pauli_term == term.pauli_term: + new_angle = merged[-1].angle + term.angle + if abs(new_angle) > atol: + merged[-1] = ExponentiatedPauliTerm(pauli_term=term.pauli_term, angle=new_angle) + else: + merged.pop() + else: + merged.append(term) + return PauliProductFormulaContainer(step_terms=merged, step_reps=1, num_qubits=self.num_qubits) + def to_json(self) -> dict[str, Any]: """Convert the PauliProductFormulaContainer to a dictionary for JSON serialization. diff --git a/python/src/qdk_chemistry/utils/qsharp/CircuitComposition.qs b/python/src/qdk_chemistry/utils/qsharp/CircuitComposition.qs new file mode 100644 index 000000000..7c21bddf1 --- /dev/null +++ b/python/src/qdk_chemistry/utils/qsharp/CircuitComposition.qs @@ -0,0 +1,54 @@ +// Copyright (c) Microsoft Corporation. All rights reserved. +// Licensed under the MIT License. See LICENSE.txt in the project root for +// license information. + +namespace QDKChemistry.Utils.CircuitComposition { + + import Std.Arrays.Subarray; + + /// Applies two operations sequentially on the same system register. + operation ApplySequential( + first : Qubit[] => Unit, + second : Qubit[] => Unit, + systems : Qubit[] + ) : Unit { + first(systems); + second(systems); + } + + /// Returns a composed operation that applies ``first`` and then ``second``. + function MakeSequentialOp(first : Qubit[] => Unit, second : Qubit[] => Unit) : Qubit[] => Unit { + ApplySequential(first, second, _) + } + + /// Returns the maximum element of the given array of integers. + function MaxInt(values : Int[]) : Int { + // Caller is responsible for not passing an empty array. + mutable max = values[0]; + for idx in 1 .. Length(values) - 1 { + let value = values[idx]; + if (value > max) { + set max = value; + } + } + return max; + } + + /// Creates a circuit for sequentially applying two operations on the same target qubits. + operation MakeSequentialCircuit( + first : Qubit[] => Unit, + second : Qubit[] => Unit, + targets : Int[] + ) : Unit { + if (Length(targets) == 0) { + // No target indices: allocate an empty register and do nothing. + use qs = Qubit[0]; + ApplySequential(first, second, qs); + } else { + // Allocate enough qubits so that all indices in 'targets' are valid. + let maxTarget = MaxInt(targets); + use qs = Qubit[1 + maxTarget]; + ApplySequential(first, second, Subarray(targets, qs)); + } + } +} diff --git a/python/src/qdk_chemistry/utils/qsharp/PauliExp.qs b/python/src/qdk_chemistry/utils/qsharp/PauliExp.qs new file mode 100644 index 000000000..c9ef598f6 --- /dev/null +++ b/python/src/qdk_chemistry/utils/qsharp/PauliExp.qs @@ -0,0 +1,92 @@ +// Copyright (c) Microsoft Corporation. All rights reserved. +// Licensed under the MIT License. See LICENSE.txt in the project root for +// license information. + +namespace QDKChemistry.Utils.PauliExp { + + import Std.Arrays.Subarray; + + /// Performs Time Evolution for a set of Pauli exponentials. + /// # Parameters + /// - `pauliExponents`: An array of arrays of Pauli operators representing the Pauli terms. + /// - `pauliCoefficients`: An array of doubles representing the coefficients for each Pauli term. + /// - `systems`: An array of qubits representing the system on which the operation acts. + /// # Returns + /// - `Unit`: The operation prepares the time evolution on the allocated qubits. + operation PauliExp( + pauliExponents : Pauli[][], + pauliCoefficients : Double[], + systems : Qubit[] + ) : Unit { + for idx in 0..Length(pauliExponents) - 1 { + let paulis = pauliExponents[idx]; + let coeff = pauliCoefficients[idx]; + Exp(paulis, -coeff, systems); + } + } + + + /// Performs repeated Time Evolution for a set of Pauli exponentials. + /// # Parameters + /// - `pauliExponents`: An array of arrays of Pauli operators representing the Pauli terms. + /// - `pauliCoefficients`: An array of doubles representing the coefficients for each Pauli term. + /// - `repetitions`: The number of times to repeat the evolution. + struct RepPauliExpParams { + pauliExponents : Pauli[][], + pauliCoefficients : Double[], + repetitions : Int, + } + + /// Performs repeated Time Evolution for a set of Pauli exponentials. + /// # Parameters + /// - `params`: A `RepPauliExpParams` struct containing the parameters for the operation. + /// - `systems`: An array of qubits representing the system on which the operation acts. + /// # Returns + /// - `Unit`: The operation prepares the repeated time evolution on the allocated qubits. + operation RepPauliExp( + params : RepPauliExpParams, + systems : Qubit[], + ) : Unit { + for i in 1..params.repetitions { + PauliExp(params.pauliExponents, params.pauliCoefficients, systems); + } + } + + /// A helper operation to create a circuit for repeated Time Evolution for a set of Pauli exponentials. + /// # Parameters + /// - `params`: A `RepPauliExpParams` struct containing the parameters for the operation. + /// - `system`: An array of integers representing the indices of the system qubits. + /// # Returns + /// - `Unit`: The operation prepares the repeated time evolution on the allocated qubits. + operation MakeRepPauliExpCircuit( + params : RepPauliExpParams, + system : Int[], + ) : Unit { + // If no system indices are provided, there is nothing to do. + if Length(system) == 0 { + return (); + } + + // Determine the maximum index in the system array to size the qubit register safely. + mutable maxIndex = system[0]; + for idx in 1..Length(system) - 1 { + let current = system[idx]; + if current > maxIndex { + set maxIndex = current; + } + } + + // Allocate enough qubits so that all indices in `system` are valid. + use qs = Qubit[maxIndex + 1]; + RepPauliExp(params, Subarray(system, qs)); + } + + /// A helper function to create a callable for repeated Time Evolution for a set of Pauli exponentials. + /// # Parameters + /// - `params`: A `RepPauliExpParams` struct containing the parameters for the operation. + /// # Returns + /// - `Qubit[] => Unit`: A callable that takes an array of system qubits, and prepares the repeated time evolution on the allocated qubits. + function MakeRepPauliExpOp(params : RepPauliExpParams) : Qubit[] => Unit { + RepPauliExp(params, _) + } +} diff --git a/python/src/qdk_chemistry/utils/qsharp/__init__.py b/python/src/qdk_chemistry/utils/qsharp/__init__.py index 3cbe37545..6a271f177 100644 --- a/python/src/qdk_chemistry/utils/qsharp/__init__.py +++ b/python/src/qdk_chemistry/utils/qsharp/__init__.py @@ -13,8 +13,10 @@ _QS_FILES = [ Path(__file__).parent / "StatePreparation.qs", + Path(__file__).parent / "CircuitComposition.qs", Path(__file__).parent / "IterativePhaseEstimation.qs", Path(__file__).parent / "ControlledPauliExp.qs", + Path(__file__).parent / "PauliExp.qs", Path(__file__).parent / "MeasurementBasis.qs", ] diff --git a/python/tests/test_evolve_and_measure.py b/python/tests/test_evolve_and_measure.py new file mode 100644 index 000000000..0c808093f --- /dev/null +++ b/python/tests/test_evolve_and_measure.py @@ -0,0 +1,123 @@ +"""Tests for EvolveAndMeasure state-preparation composition.""" + +# -------------------------------------------------------------------------------------------- +# Copyright (c) Microsoft Corporation. All rights reserved. +# Licensed under the MIT License. See LICENSE.txt in the project root for license information. +# -------------------------------------------------------------------------------------------- + +from __future__ import annotations + +from types import SimpleNamespace + +import numpy as np +import pytest + +import qdk_chemistry.algorithms.time_evolution.measure_simulation.base as measure_base +from qdk_chemistry.algorithms import create +from qdk_chemistry.algorithms.time_evolution.builder.trotter import Trotter +from qdk_chemistry.algorithms.time_evolution.circuit_mapper import PauliSequenceMapper +from qdk_chemistry.algorithms.time_evolution.measure_simulation import EvolveAndMeasure +from qdk_chemistry.data import Circuit, QubitHamiltonian + + +def test_prepend_state_prep_circuit_composes_qsharp_operations(monkeypatch: pytest.MonkeyPatch) -> None: + """The helper should compose state preparation and evolution through Q# operations.""" + algo = EvolveAndMeasure() + state_prep_op = object() + evolution_op = object() + + monkeypatch.setattr( + measure_base, + "QSHARP_UTILS", + SimpleNamespace( + CircuitComposition=SimpleNamespace( + MakeSequentialCircuit="sequential-circuit", + MakeSequentialOp=lambda first, second: ("sequential-op", first, second), + ) + ), + ) + monkeypatch.setattr( + measure_base.qsharp, + "circuit", + lambda *args: ("qsharp-circuit", args), + ) + monkeypatch.setattr( + measure_base.qsharp, + "compile", + lambda *args: ("qir", args), + ) + + state_prep = Circuit( + qir='attributes #0 = { "required_num_qubits"="1" }', + qsharp_op=state_prep_op, + encoding="jordan-wigner", + ) + evolution = Circuit( + qir='attributes #0 = { "required_num_qubits"="1" }', + qsharp_op=evolution_op, + encoding="jordan-wigner", + ) + + combined = algo._prepend_state_prep_circuit(state_prep, evolution) + + assert combined.qsharp == ( + "qsharp-circuit", + ("sequential-circuit", state_prep_op, evolution_op, [0]), + ) + assert combined.qir == ( + "qir", + ("sequential-circuit", state_prep_op, evolution_op, [0]), + ) + assert combined._qsharp_op == ("sequential-op", state_prep_op, evolution_op) + assert combined.encoding == "jordan-wigner" + + +def test_prepend_state_prep_circuit_rejects_mismatched_qubit_counts() -> None: + """The helper should reject incompatible circuit widths before composing.""" + algo = EvolveAndMeasure() + state_prep = Circuit(qir='attributes #0 = { "required_num_qubits"="1" }', qsharp_op=lambda _: None) + evolution = Circuit(qir='attributes #0 = { "required_num_qubits"="2" }', qsharp_op=lambda _: None) + + with pytest.raises(ValueError, match="same number of qubits"): + algo._prepend_state_prep_circuit(state_prep, evolution) + + +def test_prepend_state_prep_circuit_requires_qsharp_operations() -> None: + """The helper should fail fast when either circuit lacks a Q# operation handle.""" + algo = EvolveAndMeasure() + state_prep = Circuit(qasm="OPENQASM 3.0;\nqubit[1] q;\nh q[0];\n") + evolution = Circuit(qasm="OPENQASM 3.0;\nqubit[1] q;\nx q[0];\n") + + with pytest.raises(RuntimeError, match="requires Q# operations"): + algo._prepend_state_prep_circuit(state_prep, evolution) + + +def test_evolve_and_measure_eigenvalue_remains_constant() -> None: + """Run an example EvolveAndMeasure workflow.""" + hamiltonian_p = QubitHamiltonian(["ZZ"], np.array([1.0])) + hamiltonian_m = QubitHamiltonian(["ZZ"], np.array([-1.0])) + + steps = 100 + hamiltonians = [hamiltonian_p, hamiltonian_m] * (steps // 2) + time_steps = [float(t + 1) for t in range(steps)] + observable = QubitHamiltonian(["ZZ"], np.array([1.0])) + + evolution_builder = Trotter(num_divisions=1, order=1, optimize_term_ordering=True) + algo = EvolveAndMeasure() + mapper = PauliSequenceMapper() + energy_estimator = create("energy_estimator", "qdk") + circuit_executor = create("circuit_executor", "qdk_full_state_simulator") + + measurements = algo.run( + hamiltonians, + times=time_steps, + observables=[observable], + evolution_builder=evolution_builder, + circuit_mapper=mapper, + circuit_executor=circuit_executor, + energy_estimator=energy_estimator, + shots=1024, + ) + + for measurement in measurements: + assert measurement[0].energy_expectation_value == pytest.approx(1.0, abs=0.2) diff --git a/python/tests/test_time_evolution_circuit_mapper_noncontrolled.py b/python/tests/test_time_evolution_circuit_mapper_noncontrolled.py new file mode 100644 index 000000000..b65d753ed --- /dev/null +++ b/python/tests/test_time_evolution_circuit_mapper_noncontrolled.py @@ -0,0 +1,94 @@ +"""Tests for the non-controlled PauliSequenceMapper in QDK/Chemistry.""" + +# -------------------------------------------------------------------------------------------- +# Copyright (c) Microsoft Corporation. All rights reserved. +# Licensed under the MIT License. See LICENSE.txt in the project root for license information. +# -------------------------------------------------------------------------------------------- + +import json + +import numpy as np +import pytest +import qsharp +import scipy + +from qdk_chemistry.algorithms.time_evolution.circuit_mapper.pauli_sequence_mapper import PauliSequenceMapper +from qdk_chemistry.data.circuit import Circuit +from qdk_chemistry.data.time_evolution.base import TimeEvolutionUnitary +from qdk_chemistry.data.time_evolution.containers.pauli_product_formula import ( + ExponentiatedPauliTerm, + PauliProductFormulaContainer, +) +from qdk_chemistry.plugins.qiskit import QDK_CHEMISTRY_HAS_QISKIT + +from .reference_tolerances import float_comparison_absolute_tolerance, float_comparison_relative_tolerance + +if QDK_CHEMISTRY_HAS_QISKIT: + from qiskit.quantum_info import Operator + + +@pytest.fixture +def simple_unitary() -> TimeEvolutionUnitary: + """Create a simple TimeEvolutionUnitary for testing.""" + container = PauliProductFormulaContainer( + step_terms=[ + ExponentiatedPauliTerm(pauli_term={0: "X"}, angle=0.5), + ExponentiatedPauliTerm(pauli_term={1: "Z"}, angle=0.25), + ], + step_reps=2, + num_qubits=2, + ) + return TimeEvolutionUnitary(container=container) + + +class TestPauliSequenceMapperNonControlled: + """Tests for the non-controlled PauliSequenceMapper class.""" + + def test_name_and_type_name(self): + """Test mapper identity methods.""" + mapper = PauliSequenceMapper() + + assert mapper.name() == "pauli_sequence" + assert mapper.type_name() == "evolution_circuit_mapper" + + def test_run_builds_regular_unitary_circuit(self, simple_unitary): + """Test run() builds a regular (non-controlled) unitary circuit.""" + mapper = PauliSequenceMapper() + + circuit = mapper.run(simple_unitary) + + assert isinstance(circuit, Circuit) + assert isinstance(circuit.get_qsharp_circuit(), qsharp._native.Circuit) + + qsc_json = json.loads(circuit.get_qsharp_circuit().json()) + num_qubits = len(qsc_json["qubits"]) + assert num_qubits == 2 + + @pytest.mark.skipif(not QDK_CHEMISTRY_HAS_QISKIT, reason="Qiskit not available.") + def test_unitary_circuit_matrix(self, simple_unitary): + """Test that the constructed unitary circuit has the expected matrix.""" + mapper = PauliSequenceMapper() + circuit = mapper.run(simple_unitary) + + container = simple_unitary.get_container() + angle_x = container.step_terms[0].angle + angle_z = container.step_terms[1].angle + + pauli_x = np.array([[0, 1], [1, 0]], dtype=complex) + pauli_z = np.array([[1, 0], [0, -1]], dtype=complex) + identity = np.eye(2, dtype=complex) + x_0 = np.kron(identity, pauli_x) + z_1 = np.kron(pauli_z, identity) + + u_step = scipy.linalg.expm(-1j * angle_z * z_1) @ scipy.linalg.expm(-1j * angle_x * x_0) + expected_matrix = np.linalg.matrix_power(u_step, container.step_reps) + + qc = circuit.get_qiskit_circuit() + actual_matrix = Operator(qc).data + + assert np.allclose( + actual_matrix, + expected_matrix, + atol=float_comparison_absolute_tolerance, + rtol=float_comparison_relative_tolerance, + ) diff --git a/python/tests/test_time_evolution_container.py b/python/tests/test_time_evolution_container.py index 6089e7ce9..0e15d8480 100644 --- a/python/tests/test_time_evolution_container.py +++ b/python/tests/test_time_evolution_container.py @@ -128,6 +128,70 @@ def test_to_hdf5_roundtrip(self, container, tmp_path): assert restored.step_reps == container.step_reps assert len(restored.step_terms) == len(container.step_terms) + def test_combine_no_adjacent_identical(self): + """Test combine when no adjacent terms share the same Pauli string.""" + a = PauliProductFormulaContainer( + step_terms=[ + ExponentiatedPauliTerm(pauli_term={0: "X"}, angle=0.1), + ExponentiatedPauliTerm(pauli_term={1: "Z"}, angle=0.2), + ], + step_reps=2, + num_qubits=2, + ) + b = PauliProductFormulaContainer( + step_terms=[ + ExponentiatedPauliTerm(pauli_term={0: "Y"}, angle=0.3), + ExponentiatedPauliTerm(pauli_term={0: "X"}, angle=0.4), + ], + step_reps=2, + num_qubits=2, + ) + result = a.combine(b) + + # a expanded: [X, Z, X, Z], b expanded: [Y, X, Y, X] + # No adjacent duplicates anywhere, so all 8 terms survive. + assert result.step_reps == 1 + assert len(result.step_terms) == 8 + expected_angles = [0.1, 0.2, 0.1, 0.2, 0.3, 0.4, 0.3, 0.4] + for term, expected in zip(result.step_terms, expected_angles, strict=True): + assert np.isclose(term.angle, expected, atol=1e-14) + + def test_combine_with_adjacent_identical(self): + """Test combine where adjacent identical Pauli terms get merged.""" + a = PauliProductFormulaContainer( + step_terms=[ + ExponentiatedPauliTerm(pauli_term={0: "Y"}, angle=1.5), + ExponentiatedPauliTerm(pauli_term={0: "X"}, angle=0.5), + ], + step_reps=2, + num_qubits=1, + ) + b = PauliProductFormulaContainer( + step_terms=[ + ExponentiatedPauliTerm(pauli_term={0: "X"}, angle=0.7), + ExponentiatedPauliTerm(pauli_term={0: "Z"}, angle=1.5), + ], + step_reps=1, + num_qubits=1, + ) + result = a.combine(b) + + # a expanded: [Y(1.5), X(0.5), Y(1.5), X(0.5)], b expanded: [X(0.7), Z(1.5)] + # Only the two adjacent X terms at the boundary are merged into X(1.2) + assert result.step_reps == 1 + assert len(result.step_terms) == 5 + + assert result.step_terms[0].pauli_term == {0: "Y"} + assert np.isclose(result.step_terms[0].angle, 1.5, atol=1e-14) + assert result.step_terms[1].pauli_term == {0: "X"} + assert np.isclose(result.step_terms[1].angle, 0.5, atol=1e-14) + assert result.step_terms[2].pauli_term == {0: "Y"} + assert np.isclose(result.step_terms[2].angle, 1.5, atol=1e-14) + assert result.step_terms[3].pauli_term == {0: "X"} + assert np.isclose(result.step_terms[3].angle, 1.2, atol=1e-14) + assert result.step_terms[4].pauli_term == {0: "Z"} + assert np.isclose(result.step_terms[4].angle, 1.5, atol=1e-14) + def test_summary(self, container): """Test the summary generation of the container.""" summary = container.get_summary() diff --git a/python/tests/test_time_evolution_trotter.py b/python/tests/test_time_evolution_trotter.py index eed5913f8..e2446f51b 100644 --- a/python/tests/test_time_evolution_trotter.py +++ b/python/tests/test_time_evolution_trotter.py @@ -104,6 +104,30 @@ def test_multiple_trotter_steps(self): rtol=float_comparison_relative_tolerance, ) + def test_single_step_merge_identical_terms(self): + """Test construction of time evolution unitary with a single Trotter step.""" + pauli_strings = ["XII", "IXI", "XII"] + coefficients = [1.0, 1.0, 1.0] + + # Scramble the order of the terms to ensure grouping works, using a fixed seed + perm = np.random.default_rng(seed=0).permutation(len(pauli_strings)) + pauli_strings = [pauli_strings[i] for i in perm] + coefficients = [coefficients[i] for i in perm] + + hamiltonian = QubitHamiltonian(pauli_strings=pauli_strings, coefficients=coefficients) + builder = Trotter(num_divisions=1, optimize_term_ordering=True) + unitary = builder.run(hamiltonian, time=1) + + assert isinstance(unitary, TimeEvolutionUnitary) + container = unitary.get_container() + + assert isinstance(container, PauliProductFormulaContainer) + assert container.num_qubits == 3 + assert container.step_reps == 1 + # After merging identical Pauli strings, 2 unique terms remain (XII with + # coeff=2 and IXI with coeff=1). Raw Pauli terms are emitted directly. + assert len(container.step_terms) == 2 + def test_basic_decomposition(self): """Test basic decomposition of a qubit Hamiltonian.""" builder = Trotter() @@ -834,3 +858,53 @@ def test_angle_scaling_with_auto_steps_higher_order(self): atol=float_comparison_absolute_tolerance, rtol=float_comparison_relative_tolerance, ) + + +class TestOptimizeTermOrdering: + """Tests for the optimize_term_ordering method.""" + + def test_optimize_term_ordering_does_nothing_when_false(self): + """Test that optimize_term_ordering does nothing when optimize_term_ordering is False.""" + pauli_strings = ["ZZII", "XIII", "IZZI", "IXII", "IIZZ", "IIXI", "ZIIZ", "IIIX"] + coefficients = [1.0] * 8 + hamiltonian = QubitHamiltonian( + pauli_strings=pauli_strings, + coefficients=coefficients, + ) + builder = Trotter(num_divisions=1, order=1, optimize_term_ordering=False) + t = 1.0 + terms = builder.run(hamiltonian, time=t).get_container().step_terms + + for idx, term in enumerate(terms): + assert term.pauli_term == builder._pauli_label_to_map(pauli_strings[idx]) + assert term.angle == coefficients[idx] * t + + def test_optimize_term_ordering_groups_when_true(self): + """Test that optimize_term_ordering groups commuting terms into parallelizable layers.""" + pauli_strings = ["ZZII", "XIII", "IZZI", "IXII", "IIZZ", "IIXI", "ZIIZ", "IIIX"] + coefficients = [1.0] * 8 + hamiltonian = QubitHamiltonian( + pauli_strings=pauli_strings, + coefficients=coefficients, + ) + builder = Trotter(num_divisions=1, order=1, optimize_term_ordering=True) + t = 1.0 + terms = builder.run(hamiltonian, time=t).get_container().step_terms + + # Convert terms back to label strings for grouping + def term_to_label(term): + num_qubits = hamiltonian.num_qubits + chars = ["I"] * num_qubits + for qubit_idx, pauli_op in term.pauli_term.items(): + chars[num_qubits - 1 - qubit_idx] = pauli_op + return "".join(chars) + + term_labels = [term_to_label(t) for t in terms] + + # Raw Pauli terms are emitted directly (no Clifford sandwich). + # ZZ terms stay as their original labels; X terms stay as single-qubit X. + pauli_zz_labels = {"ZZII", "IIZZ", "IZZI", "ZIIZ"} + pauli_x_labels = {"XIII", "IXII", "IIXI", "IIIX"} + + assert len(terms) == len(pauli_zz_labels) + len(pauli_x_labels) + assert set(term_labels) == pauli_zz_labels | pauli_x_labels