diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index daa169c3..e0769ec1 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -180,8 +180,8 @@ jobs: - name: Install and Test wheel with extras run: | pip install setuptools - bash bin/install_wheel_extras.sh dist --type wheel --extra test --extra cli --extra pulse - python -m pytest --cov=pyqasm tests/ --cov-report=html --cov-report=xml --cov-report=term + bash bin/install_wheel_extras.sh dist --type wheel --extra test --extra test-sim --extra cli --extra pulse + python -m pytest --cov=pyqasm tests/ -m "not benchmark" --cov-report=html:build/coverage/html --cov-report=xml:build/coverage/coverage.xml --cov-report=term - name: Upload coverage to Codecov if: ${{ matrix.python == 311 }} @@ -190,4 +190,4 @@ jobs: token: ${{ secrets.CODECOV_TOKEN }} fail_ci_if_error: false files: ./build/coverage/coverage.xml - verbose: true \ No newline at end of file + verbose: true diff --git a/.github/workflows/pre-release.yml b/.github/workflows/pre-release.yml index c3b99e11..b111d55b 100644 --- a/.github/workflows/pre-release.yml +++ b/.github/workflows/pre-release.yml @@ -176,4 +176,4 @@ jobs: uses: pypa/gh-action-pypi-publish@release/v1 with: user: __token__ - password: ${{ secrets.PYPI_API_TOKEN }} \ No newline at end of file + password: ${{ secrets.PYPI_API_TOKEN }} diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index d660a47e..e244cef3 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -167,4 +167,4 @@ jobs: uses: pypa/gh-action-pypi-publish@release/v1 with: user: __token__ - password: ${{ secrets.PYPI_API_TOKEN }} \ No newline at end of file + password: ${{ secrets.PYPI_API_TOKEN }} diff --git a/.github/workflows/test-release.yml b/.github/workflows/test-release.yml index 838f1aa9..0bedb1f7 100644 --- a/.github/workflows/test-release.yml +++ b/.github/workflows/test-release.yml @@ -166,4 +166,4 @@ jobs: with: repository-url: https://test.pypi.org/legacy/ user: __token__ - password: ${{ secrets.TEST_PYPI_API_TOKEN }} \ No newline at end of file + password: ${{ secrets.TEST_PYPI_API_TOKEN }} diff --git a/.gitignore b/.gitignore index 09e6a753..c9961e83 100644 --- a/.gitignore +++ b/.gitignore @@ -8,6 +8,7 @@ __pycache__/ # Cython build files src/pyqasm/accelerate/linalg.c +src/pyqasm/accelerate/sv_sim.c # Distribution / packaging .Python diff --git a/CHANGELOG.md b/CHANGELOG.md index 5dde4024..e7963eb3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,7 @@ Types of changes: ## Unreleased ### Added +- Added a statevector simulator as a new `pyqasm.simulator` subpackage, backed by a Cython/OpenMP-capable kernel (`pyqasm.accelerate.sv_sim`). `Simulator.run()` accepts an OpenQASM 3 string or an unrolled `QasmModule` and returns a `SimulatorResult` with the final statevector, probabilities, and sampled measurement counts. Ships with cross-validation tests against qiskit (`test-sim` extra), an optional `simulation` extra (numba-accelerated preprocessing helpers), and a `benchmarks/` suite. ([#316](https://github.com/qBraid/pyqasm/pull/316)) - Added support for the `c3x` (3-controlled X) and `rc3x`/`rcccx` (relative-phase 3-controlled X) gates, decomposed into basis gates following qiskit's `C3XGate`/`RC3XGate` definitions. Also extended the `ctrl @` modifier chain so that 3- and 4-control stacks on `x` (e.g. `ctrl @ ctrl @ ctrl @ x`, `ctrl(4) @ x`) resolve to `c3x`/`c4x`. ([#320](https://github.com/qBraid/pyqasm/pull/320)) ### Improved / Modified @@ -27,6 +28,7 @@ Types of changes: ### Fixed - Fixed `reset` on a physical qubit rewriting the operand to the internal pulse register, e.g. `reset $2;` unrolled to `reset __PYQASM_QUBITS__[2];`. That names a register the program never declares, so the unrolled output did not round-trip through `dumps()`/`loads()`, and the qubit was never registered (a program whose only operation was `reset $3;` reported `num_qubits == 0`). Physical qubits are now kept as-is in plain QASM programs, matching how gate and measurement operands already treat them; the rename is still applied for OpenPulse programs, where the pulse visitor expects it. ([#325](https://github.com/qBraid/pyqasm/pull/325)) - Fixed `measure` and `reset` on a user register whose name merely starts with the reserved internal register name (e.g. `__PYQASM_QUBITS__foo`) being mistaken for the internal register itself. Such statements were short-circuited out of unrolling and emitted verbatim, so `c = measure __PYQASM_QUBITS__foo;` was never expanded per qubit and `reset __PYQASM_QUBITS__foo;` silently reset nothing. The register is now matched on its exact name, or on the name followed by an index or slice. ([#325](https://github.com/qBraid/pyqasm/pull/325)) +- Fixed the statevector simulator's `crz` fast path, which implemented a controlled-phase gate (`diag(1, 1, 1, e^{i*theta/2})`) instead of a true controlled-Rz (`diag(1, 1, e^{-i*theta/2}, e^{i*theta/2})`); it dropped the phase on the `|control=1, target=0>` amplitude. It now applies the full `Rz` to the target via the controlled-gate kernel. ([#316](https://github.com/qBraid/pyqasm/pull/316)) - Fixed the `c4x` (4-controlled X) gate, which previously raised `TypeError: c4x_gate() takes 4 positional arguments but 5 were given` because it was declared with four parameters for a five-qubit gate. It is now implemented via qiskit's structured `rc3x`/`c3sx`/`cphaseshift` decomposition. ([#320](https://github.com/qBraid/pyqasm/pull/320)) - Added `inv @` (inverse modifier) support for the multi-controlled-X family: `inv @ c3x` / `inv @ c4x` resolve to the (self-inverse) gate, and `inv @ rc3x` / `inv @ rcccx` resolve to the correct relative-phase dagger. These previously raised `Unsupported / undeclared QASM operation`. ([#320](https://github.com/qBraid/pyqasm/pull/320)) - Fixed the `ctrl @` modifier not resolving gate aliases: `ctrl @ toffoli` / `ctrl @ ccnot` (aliases of `ccx`) and `ctrl @ cnot` / `ctrl @ CX` (aliases of `cx`) now escalate controls identically to their canonical gate instead of raising `Unsupported controlled QASM operation`. ([#320](https://github.com/qBraid/pyqasm/pull/320)) diff --git a/MANIFEST.in b/MANIFEST.in index 5fcba451..106cccb1 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -2,4 +2,8 @@ include README.md recursive-include src/pyqasm *.py *.pyx include src/pyqasm/py.typed -exclude tests/* \ No newline at end of file +# Never ship Cython-generated C: sdists compile from .pyx (Cython is a build +# requirement), and stale generated C must not mask the .pyx sources. +global-exclude *.c + +exclude tests/* diff --git a/benchmarks/README.md b/benchmarks/README.md new file mode 100644 index 00000000..9ec43945 --- /dev/null +++ b/benchmarks/README.md @@ -0,0 +1,78 @@ +# Benchmarks + +Scripts and reference results for the PyQASM statevector simulator +(`pyqasm.simulator.Simulator`). + +## Contents + +| File | What it is | +| --- | --- | +| `bench_simulator.py` | Benchmarks PyQASM against Qiskit Aer, Cirq, and PennyLane Lightning on shared random and QFT circuits, with cross-simulator correctness checks. | +| `plot_benchmarks.py` | Runs the same benchmarks and renders log-scale comparison plots (`bench_random.png`, `bench_qft.png`). | +| `profile_simulator.py` | Profiles the PyQASM simulator itself: cProfile, section timing (preprocess vs simulation vs post-processing), or per-kernel timing. | +| `bench_random.png`, `bench_qft.png` | Committed comparison plots produced by `plot_benchmarks.py`. | +| `bench_evolved.json`, `bench_qft_evolved.png`, `bench_random_evolved.png` | Historical before/after results for the evolved preprocessing change (see provenance note below). | + +## Requirements + +All scripts need only `pyqasm` (installed from this repo) for the PyQASM +paths. Additional dependencies: + +- `bench_simulator.py`: optionally `qiskit` + `qiskit-aer`, `cirq`, and + `pennylane` (Lightning) for the comparison columns, and `tabulate` for the + summary table. Every one of these is optional — missing simulators are + reported as `skipped: not installed` and the script runs the + remaining backends (worst case, PyQASM only). +- `plot_benchmarks.py`: `matplotlib`, plus the same optional simulators as + `bench_simulator.py`. Plots include only the installed backends. +- `profile_simulator.py`: no extra dependencies. + +None of these are declared as a package extra; install them ad hoc, e.g. +`pip install qiskit qiskit-aer cirq pennylane tabulate matplotlib`. + +## Running + +From the repository root: + +```bash +# Full benchmark table (up to 22 qubits; minutes of runtime) +python benchmarks/bench_simulator.py + +# Smoke test: small circuits, 2 repeats +python benchmarks/bench_simulator.py --quick + +# Regenerate bench_random.png / bench_qft.png in benchmarks/ +python benchmarks/plot_benchmarks.py + +# Quick plot run, writing PNGs somewhere else +python benchmarks/plot_benchmarks.py --quick --out-dir /tmp/bench-plots + +# Profile the simulator (modes: section | cprofile | kernel) +python benchmarks/profile_simulator.py --mode section --qubits 16 --depth 200 +python benchmarks/profile_simulator.py --mode cprofile --qubits 12 --depth 100 +python benchmarks/profile_simulator.py --mode kernel --qubits 12 --depth 100 +``` + +Timing covers simulation only; circuit construction/transpilation is done in +each backend's `prepare_*` step and excluded from the measured time. Reported +numbers are the median of `N_REPEATS` runs (default 5; `--quick` uses 2, +`--repeats N` overrides). + +## Provenance of the `*_evolved` files + +`bench_evolved.json`, `bench_qft_evolved.png`, and `bench_random_evolved.png` +were produced during development of the evolved no-cache preprocessing change +(PR #316). They compare two states of this branch: + +- **baseline**: the branch state *before* the evolved-preprocessing commit + (the previous statevector preprocessing implementation); +- **candidate**: the branch with the evolved no-cache preprocessing. + +They are retained for historical reference only. They are **not** regenerated +by any script in this tree: the baseline implementation no longer exists in +the working tree, so the before/after comparison cannot be reproduced from +the current sources alone. To reproduce something comparable, check out the +commit preceding the evolved-preprocessing commit and benchmark both +revisions with `bench_simulator.py`. Current-code performance can always be +re-measured with `bench_simulator.py` / `plot_benchmarks.py`, which +regenerate `bench_random.png` and `bench_qft.png`. diff --git a/benchmarks/bench_evolved.json b/benchmarks/bench_evolved.json new file mode 100644 index 00000000..941d4b03 --- /dev/null +++ b/benchmarks/bench_evolved.json @@ -0,0 +1,183 @@ +{ + "benchmark": "statevector preprocessing no-cache comparison", + "baseline": "PyQASM sv branch before evolved preprocessing commit", + "candidate": "PyQASM sv branch with evolved no-cache preprocessing", + "n_repeats": 5, + "pyqasm_num_threads": 1, + "notes": [ + "No QASM string memoization is included in the candidate.", + "Both series use pre-unrolled QasmModule inputs, matching the branch benchmark convention." + ], + "series": { + "random": [ + { + "num_qubits": 4, + "depth": 20, + "pyqasm_sv_ms": 0.06092991679906845, + "evolved_no_cache_ms": 0.06546999793499708, + "speedup": 0.9306540204807044, + "correct": true + }, + { + "num_qubits": 6, + "depth": 40, + "pyqasm_sv_ms": 0.1580890966579318, + "evolved_no_cache_ms": 0.1151899341493845, + "speedup": 1.372421104546456, + "correct": true + }, + { + "num_qubits": 8, + "depth": 60, + "pyqasm_sv_ms": 0.16204908024519682, + "evolved_no_cache_ms": 0.14194007962942123, + "speedup": 1.1416724625509327, + "correct": true + }, + { + "num_qubits": 10, + "depth": 80, + "pyqasm_sv_ms": 0.24963996838778257, + "evolved_no_cache_ms": 0.21547905635088682, + "speedup": 1.1585347208002805, + "correct": true + }, + { + "num_qubits": 12, + "depth": 100, + "pyqasm_sv_ms": 0.458328053355217, + "evolved_no_cache_ms": 0.41990901809185743, + "speedup": 1.0914937131808757, + "correct": true + }, + { + "num_qubits": 14, + "depth": 150, + "pyqasm_sv_ms": 1.5235739992931485, + "evolved_no_cache_ms": 1.3638049131259322, + "speedup": 1.117149516495739, + "correct": true + }, + { + "num_qubits": 16, + "depth": 200, + "pyqasm_sv_ms": 7.1852849796414375, + "evolved_no_cache_ms": 5.937180016189814, + "speedup": 1.210218480835721, + "correct": true + }, + { + "num_qubits": 18, + "depth": 250, + "pyqasm_sv_ms": 35.15114798210561, + "evolved_no_cache_ms": 28.548280941322446, + "speedup": 1.231287728124666, + "correct": true + }, + { + "num_qubits": 20, + "depth": 300, + "pyqasm_sv_ms": 162.53185598179698, + "evolved_no_cache_ms": 126.72401999589056, + "speedup": 1.282565499319447, + "correct": true + }, + { + "num_qubits": 22, + "depth": 400, + "pyqasm_sv_ms": 1289.0826460206881, + "evolved_no_cache_ms": 1231.460615992546, + "speedup": 1.0467916141854843, + "correct": true + } + ], + "qft": [ + { + "num_qubits": 4, + "depth": "QFT", + "pyqasm_sv_ms": 0.20883907563984394, + "evolved_no_cache_ms": 0.16430008690804243, + "speedup": 1.2710831720784765, + "correct": true + }, + { + "num_qubits": 6, + "depth": "QFT", + "pyqasm_sv_ms": 0.4983479157090187, + "evolved_no_cache_ms": 0.315018929541111, + "speedup": 1.5819618091997316, + "correct": true + }, + { + "num_qubits": 8, + "depth": "QFT", + "pyqasm_sv_ms": 0.9089269442483783, + "evolved_no_cache_ms": 0.5902479169890285, + "speedup": 1.5399070764789728, + "correct": true + }, + { + "num_qubits": 10, + "depth": "QFT", + "pyqasm_sv_ms": 1.5004950109869242, + "evolved_no_cache_ms": 1.0477059986442327, + "speedup": 1.4321718238977499, + "correct": true + }, + { + "num_qubits": 12, + "depth": "QFT", + "pyqasm_sv_ms": 2.5123610394075513, + "evolved_no_cache_ms": 1.782613922841847, + "speedup": 1.409369133279482, + "correct": true + }, + { + "num_qubits": 14, + "depth": "QFT", + "pyqasm_sv_ms": 5.448751035146415, + "evolved_no_cache_ms": 4.14173596072942, + "speedup": 1.3155718005226993, + "correct": true + }, + { + "num_qubits": 16, + "depth": "QFT", + "pyqasm_sv_ms": 18.029507948085666, + "evolved_no_cache_ms": 15.524596092291176, + "speedup": 1.1613511772482323, + "correct": true + }, + { + "num_qubits": 18, + "depth": "QFT", + "pyqasm_sv_ms": 77.3717820411548, + "evolved_no_cache_ms": 68.92851099837571, + "speedup": 1.1224931587885028, + "correct": true + }, + { + "num_qubits": 20, + "depth": "QFT", + "pyqasm_sv_ms": 382.6392919290811, + "evolved_no_cache_ms": 349.38782698009163, + "speedup": 1.0951706452866319, + "correct": true + } + ] + }, + "summary": { + "random": { + "geomean_speedup": 1.1521318658620598, + "min_speedup": 0.9306540204807044, + "max_speedup": 1.372421104546456, + "all_correct": true + }, + "qft": { + "geomean_speedup": 1.3147632402156637, + "min_speedup": 1.0951706452866319, + "max_speedup": 1.5819618091997316, + "all_correct": true + } + } +} diff --git a/benchmarks/bench_qft.png b/benchmarks/bench_qft.png new file mode 100644 index 00000000..78d2eb64 Binary files /dev/null and b/benchmarks/bench_qft.png differ diff --git a/benchmarks/bench_qft_evolved.png b/benchmarks/bench_qft_evolved.png new file mode 100644 index 00000000..2dab0d29 Binary files /dev/null and b/benchmarks/bench_qft_evolved.png differ diff --git a/benchmarks/bench_random.png b/benchmarks/bench_random.png new file mode 100644 index 00000000..a783fa11 Binary files /dev/null and b/benchmarks/bench_random.png differ diff --git a/benchmarks/bench_random_evolved.png b/benchmarks/bench_random_evolved.png new file mode 100644 index 00000000..27228cc4 Binary files /dev/null and b/benchmarks/bench_random_evolved.png differ diff --git a/benchmarks/bench_simulator.py b/benchmarks/bench_simulator.py new file mode 100644 index 00000000..45073828 --- /dev/null +++ b/benchmarks/bench_simulator.py @@ -0,0 +1,516 @@ +#!/usr/bin/env python3 +# Copyright 2025 qBraid +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Benchmark: PyQASM statevector simulator vs Qiskit Aer, Cirq, +and PennyLane Lightning. + +Each simulator builds circuits programmatically from a shared gate list so the +comparison is apples-to-apples. Timing covers only simulation (circuit +construction is excluded). + +Third-party simulators (qiskit + qiskit-aer, cirq, pennylane) are optional: +any that are not installed are skipped with a message, so the script runs +with only pyqasm installed. + +Usage: + python benchmarks/bench_simulator.py [--quick] [--repeats N] +""" + +import argparse +import importlib.util +import math +import random +import time +from collections.abc import Callable +from typing import Any + +import numpy as np + +try: + from tabulate import tabulate + + HAS_TABULATE = True +except ImportError: + HAS_TABULATE = False + + +def _has_module(name: str) -> bool: + """Return whether the named module is importable without importing it.""" + return importlib.util.find_spec(name) is not None + + +HAS_QISKIT = _has_module("qiskit") and _has_module("qiskit_aer") +HAS_CIRQ = _has_module("cirq") +HAS_PENNYLANE = _has_module("pennylane") + +# --------------------------------------------------------------------------- +# Shared circuit representation +# --------------------------------------------------------------------------- + +# Gate spec: (name, qubits_tuple, params_tuple) +# name: "h", "x", "y", "z", "s", "t", "rx", "ry", "rz", "cx", "cy", "cz", +# "swap", "crz" +# qubits_tuple: (target,) for 1q gates, (control, target) for 2q gates +# params_tuple: () or (angle,) + +SINGLE_GATE_NAMES = ["h", "x", "y", "z", "s", "t", "rx", "ry", "rz"] +TWO_GATE_NAMES = ["cx", "cy", "cz", "swap"] + +# (name, qubit_indices, params) — see the comment block above. +GateSpec = tuple[str, tuple[int, ...], tuple[float, ...]] + + +def generate_random_gates( + num_qubits: int, depth: int, seed: int = 42 +) -> tuple[int, list[GateSpec]]: + """Return ``(num_qubits, gate_list)`` for a seeded random circuit of ``depth`` gates.""" + rng = random.Random(seed) + gates: list[GateSpec] = [] + for _ in range(depth): + if num_qubits >= 2 and rng.random() < 0.4: + name = rng.choice(TWO_GATE_NAMES) + q0, q1 = rng.sample(range(num_qubits), 2) + gates.append((name, (q0, q1), ())) + else: + name = rng.choice(SINGLE_GATE_NAMES) + q = rng.randint(0, num_qubits - 1) + if name == "rx": + gates.append((name, (q,), (1.0,))) + elif name == "ry": + gates.append((name, (q,), (0.5,))) + elif name == "rz": + gates.append((name, (q,), (0.3,))) + else: + gates.append((name, (q,), ())) + return num_qubits, gates + + +def generate_qft_gates(num_qubits: int) -> tuple[int, list[GateSpec]]: + """Return ``(num_qubits, gate_list)`` for a QFT circuit on ``num_qubits`` qubits.""" + gates: list[GateSpec] = [] + for i in range(num_qubits): + gates.append(("h", (i,), ())) + for j in range(i + 1, num_qubits): + k = j - i + angle = math.pi / (2**k) + gates.append(("crz", (j, i), (angle,))) + for i in range(num_qubits // 2): + gates.append(("swap", (i, num_qubits - 1 - i), ())) + return num_qubits, gates + + +# --------------------------------------------------------------------------- +# Gate-list → QASM 3 string (for PyQASM / Qiskit) +# --------------------------------------------------------------------------- + + +def gates_to_qasm(num_qubits: int, gates: list[GateSpec]) -> str: + """Render a shared gate list as an OpenQASM 3 program string.""" + lines = [ + "OPENQASM 3;", + 'include "stdgates.inc";', + f"qubit[{num_qubits}] q;", + ] + for name, qubits, params in gates: + if len(qubits) == 1: + q = qubits[0] + if params: + lines.append(f"{name}({params[0]}) q[{q}];") + else: + lines.append(f"{name} q[{q}];") + else: + q0, q1 = qubits + if params: + lines.append(f"{name}({params[0]}) q[{q0}], q[{q1}];") + else: + lines.append(f"{name} q[{q0}], q[{q1}];") + return "\n".join(lines) + + +# --------------------------------------------------------------------------- +# Simulator wrappers +# --------------------------------------------------------------------------- + +N_REPEATS = 5 + + +def _median_time(fn: Callable[[], Any], n_repeats: int | None = None) -> tuple[float, Any]: + """Run ``fn()`` ``n_repeats`` times (default ``N_REPEATS``) and return + ``(median_seconds, last_result)``.""" + if n_repeats is None: + n_repeats = N_REPEATS + times = [] + result = None + for _ in range(n_repeats): + start = time.perf_counter() + result = fn() + times.append(time.perf_counter() - start) + return float(np.median(times)), result + + +def _align_global_phase(reference: np.ndarray, candidate: np.ndarray) -> np.ndarray: + """Phase-align candidate statevector to reference (global phase is unphysical).""" + idx = int(np.argmax(np.abs(reference))) + if np.abs(reference[idx]) < 1e-12 or np.abs(candidate[idx]) < 1e-12: + return candidate + phase = candidate[idx] / reference[idx] + return candidate / (phase / np.abs(phase)) + + +# -- PyQASM ------------------------------------------------------------------ + + +def prepare_pyqasm(num_qubits: int, gates: list[GateSpec]) -> tuple[Any, Any]: + """Load and unroll the circuit; return ``(module, Simulator)`` ready to run.""" + from pyqasm import loads as pyqasm_loads + from pyqasm.simulator import Simulator + + qasm = gates_to_qasm(num_qubits, gates) + module = pyqasm_loads(qasm) + module.unroll() + module.remove_idle_qubits() + sim = Simulator(seed=0) + return module, sim + + +def bench_pyqasm(module: Any, sim: Any) -> tuple[float, Any]: + """Time PyQASM simulation of a prepared module; return ``(median_s, statevector)``.""" + + def run(): + return sim.run(module, shots=0).final_statevector + + return _median_time(run) + + +# -- Qiskit Aer --------------------------------------------------------------- + + +def prepare_qiskit(num_qubits: int, gates: list[GateSpec]) -> tuple[Any, Any]: + """Build and transpile the circuit; return ``(compiled_circuit, AerSimulator)``.""" + from qiskit import transpile + from qiskit.qasm3 import loads as qiskit_loads + from qiskit_aer import AerSimulator + + qasm = gates_to_qasm(num_qubits, gates) + backend = AerSimulator(method="statevector") + circuit = qiskit_loads(qasm) + circuit.save_statevector() + compiled = transpile(circuit, backend, optimization_level=0) + return compiled, backend + + +def bench_qiskit(compiled: Any, backend: Any) -> tuple[float, Any]: + """Time Qiskit Aer simulation; return ``(median_s, statevector)``.""" + + def run(): + job = backend.run(compiled) + result = job.result() + return np.asarray(result.get_statevector(compiled)) + + return _median_time(run) + + +# -- Cirq ---------------------------------------------------------------------- + + +def _build_cirq_circuit(num_qubits: int, gates: list[GateSpec]) -> tuple[Any, Any]: + """Translate the gate list to cirq; return ``(Circuit, line_qubits)``.""" + import cirq + + qubits = cirq.LineQubit.range(num_qubits) + ops = [] + for name, qubit_indices, params in gates: + if name == "h": + ops.append(cirq.H(qubits[qubit_indices[0]])) + elif name == "x": + ops.append(cirq.X(qubits[qubit_indices[0]])) + elif name == "y": + ops.append(cirq.Y(qubits[qubit_indices[0]])) + elif name == "z": + ops.append(cirq.Z(qubits[qubit_indices[0]])) + elif name == "s": + ops.append(cirq.S(qubits[qubit_indices[0]])) + elif name == "t": + ops.append(cirq.T(qubits[qubit_indices[0]])) + elif name == "rx": + ops.append(cirq.rx(params[0])(qubits[qubit_indices[0]])) + elif name == "ry": + ops.append(cirq.ry(params[0])(qubits[qubit_indices[0]])) + elif name == "rz": + ops.append(cirq.rz(params[0])(qubits[qubit_indices[0]])) + elif name == "cx": + ops.append(cirq.CNOT(qubits[qubit_indices[0]], qubits[qubit_indices[1]])) + elif name == "cy": + ops.append( + cirq.ControlledGate(cirq.Y).on(qubits[qubit_indices[0]], qubits[qubit_indices[1]]) + ) + elif name == "cz": + ops.append(cirq.CZ(qubits[qubit_indices[0]], qubits[qubit_indices[1]])) + elif name == "swap": + ops.append(cirq.SWAP(qubits[qubit_indices[0]], qubits[qubit_indices[1]])) + elif name == "crz": + ops.append( + cirq.ControlledGate(cirq.rz(params[0])).on( + qubits[qubit_indices[0]], qubits[qubit_indices[1]] + ) + ) + return cirq.Circuit(ops), qubits + + +def prepare_cirq(num_qubits: int, gates: list[GateSpec]) -> tuple[Any, Any, Any]: + """Build the cirq circuit; return ``(circuit, Simulator, qubits)``.""" + import cirq + + circuit, qubits = _build_cirq_circuit(num_qubits, gates) + sim = cirq.Simulator(dtype=np.complex128) + return circuit, sim, qubits + + +def bench_cirq(circuit: Any, sim: Any, qubits: Any) -> tuple[float, Any]: + """Time cirq simulation; return ``(median_s, statevector)``.""" + + def run(): + result = sim.simulate(circuit, qubit_order=qubits) + return result.final_state_vector + + return _median_time(run) + + +# -- PennyLane Lightning ------------------------------------------------------- + + +def _build_pennylane_fn(num_qubits: int, gates: list[GateSpec]) -> Callable[[], Any]: + """Translate the gate list to a PennyLane Lightning QNode returning the state.""" + import pennylane as qml + + dev = qml.device("lightning.qubit", wires=num_qubits) + + gate_list = list(gates) # capture for closure + + @qml.qnode(dev) + def circuit(): + for name, qubit_indices, params in gate_list: + if name == "h": + qml.Hadamard(qubit_indices[0]) + elif name == "x": + qml.PauliX(qubit_indices[0]) + elif name == "y": + qml.PauliY(qubit_indices[0]) + elif name == "z": + qml.PauliZ(qubit_indices[0]) + elif name == "s": + qml.S(qubit_indices[0]) + elif name == "t": + qml.T(qubit_indices[0]) + elif name == "rx": + qml.RX(params[0], qubit_indices[0]) + elif name == "ry": + qml.RY(params[0], qubit_indices[0]) + elif name == "rz": + qml.RZ(params[0], qubit_indices[0]) + elif name == "cx": + qml.CNOT([qubit_indices[0], qubit_indices[1]]) + elif name == "cy": + qml.CY([qubit_indices[0], qubit_indices[1]]) + elif name == "cz": + qml.CZ([qubit_indices[0], qubit_indices[1]]) + elif name == "swap": + qml.SWAP([qubit_indices[0], qubit_indices[1]]) + elif name == "crz": + qml.CRZ(params[0], [qubit_indices[0], qubit_indices[1]]) + return qml.state() + + return circuit + + +def prepare_pennylane(num_qubits: int, gates: list[GateSpec]) -> tuple[Callable[[], Any]]: + """Build the Lightning QNode; return a 1-tuple ``(circuit_fn,)``.""" + circuit_fn = _build_pennylane_fn(num_qubits, gates) + return (circuit_fn,) + + +def bench_pennylane(circuit_fn: Callable[[], Any]) -> tuple[float, Any]: + """Time PennyLane Lightning simulation; return ``(median_s, statevector)``.""" + + def run(): + return np.asarray(circuit_fn()) + + return _median_time(run) + + +# --------------------------------------------------------------------------- +# Simulator registry (only installed backends) +# --------------------------------------------------------------------------- + +# name -> (prepare_fn, bench_fn, statevector_is_big_endian) +_ALL_SIMULATORS: list[tuple[str, bool, tuple[Callable, Callable, bool]]] = [ + ("PyQASM", True, (prepare_pyqasm, bench_pyqasm, False)), + ("Qiskit Aer", HAS_QISKIT, (prepare_qiskit, bench_qiskit, False)), + ("Cirq", HAS_CIRQ, (prepare_cirq, bench_cirq, True)), + ("PennyLane Lightning", HAS_PENNYLANE, (prepare_pennylane, bench_pennylane, True)), +] + + +def available_simulators(verbose: bool = False) -> dict[str, tuple[Callable, Callable, bool]]: + """Return ``{name: (prepare_fn, bench_fn, big_endian)}`` for installed backends, + optionally printing a skip message for each missing one.""" + sims = {} + for name, installed, entry in _ALL_SIMULATORS: + if installed: + sims[name] = entry + elif verbose: + print(f"skipped: {name} not installed") + return sims + + +# --------------------------------------------------------------------------- +# Qubit ordering helpers +# --------------------------------------------------------------------------- + + +def _reverse_endian(sv: np.ndarray, num_qubits: int) -> np.ndarray: + """Convert big-endian statevector to little-endian (or vice versa).""" + return sv.reshape([2] * num_qubits).transpose(range(num_qubits - 1, -1, -1)).ravel() + + +# --------------------------------------------------------------------------- +# Main +# --------------------------------------------------------------------------- + + +def run_benchmarks( + circuit_name: str, + configs: list[tuple], + generator_fn: Callable[..., tuple[int, list[GateSpec]]], + simulators: dict[str, tuple[Callable, Callable, bool]], +) -> list[list]: + """Benchmark every available simulator on each config, print a table of + median timings with a phase-invariant correctness check against PyQASM, + and return the table rows.""" + headers = ["Qubits", "Depth"] + [f"{name} (ms)" for name in simulators] + ["Correct"] + rows = [] + + print(f"\n{'='*90}") + print(f" {circuit_name} circuits (median of {N_REPEATS} runs)") + print(f"{'='*90}") + + for config in configs: + if len(config) == 2: + num_qubits, depth = config + nq, gates = generator_fn(num_qubits, depth) + else: + num_qubits = config[0] + depth = None + nq, gates = generator_fn(num_qubits) + + depth_str = str(depth) if depth is not None else "QFT" + + # --- Prepare and benchmark each available simulator --- + times = {} + statevectors = {} + for name, (prepare_fn, bench_fn, big_endian) in simulators.items(): + args = prepare_fn(nq, gates) + t, sv = bench_fn(*args) + times[name] = t + statevectors[name] = _reverse_endian(sv, nq) if big_endian else sv + + # --- Correctness: compare everything against PyQASM (little-endian) --- + # Compare up to global phase: simulators may legitimately differ by an + # overall e^{i*phi}, which is physically irrelevant. + sv_pyqasm = statevectors["PyQASM"] + checks = { + name: np.allclose(sv_pyqasm, _align_global_phase(sv_pyqasm, sv), atol=1e-10) + for name, sv in statevectors.items() + if name != "PyQASM" + } + if not checks: + status = "n/a (PyQASM only)" + elif all(checks.values()): + status = "PASS" + else: + status = "FAIL " + ",".join(k for k, v in checks.items() if not v) + + rows.append([nq, depth_str] + [f"{times[name]*1000:.2f}" for name in simulators] + [status]) + + timings = " ".join(f"{name}={times[name]*1000:8.2f}" for name in simulators) + print(f" n={nq:2d} depth={depth_str:>4s} {timings} ms {status}") + + print() + if HAS_TABULATE: + print(tabulate(rows, headers=headers, tablefmt="github")) + else: + print(" ".join(headers)) + for row in rows: + print(" ".join(str(cell) for cell in row)) + return rows + + +def _positive_int(value: str) -> int: + """argparse type: accept only strictly positive integers.""" + parsed = int(value) + if parsed <= 0: + raise argparse.ArgumentTypeError(f"must be a positive integer, got {value!r}") + return parsed + + +def main() -> None: + """Parse CLI options and run the random-circuit and QFT benchmark suites.""" + global N_REPEATS + + parser = argparse.ArgumentParser(description="Benchmark PyQASM against other simulators") + parser.add_argument( + "--quick", + action="store_true", + help="Small qubit counts and fewer repeats (smoke test).", + ) + parser.add_argument( + "--repeats", + type=_positive_int, + default=None, + help=f"Number of timed repeats per config (default {N_REPEATS}, or 2 with --quick).", + ) + args = parser.parse_args() + + if args.repeats is not None: + N_REPEATS = args.repeats + elif args.quick: + N_REPEATS = 2 + + if args.quick: + random_configs = [(2, 10), (4, 20), (6, 30)] + qft_configs = [(2,), (4,), (6,)] + else: + random_configs = [ + (4, 20), + (6, 40), + (8, 60), + (10, 80), + (12, 100), + (14, 150), + (16, 200), + (18, 250), + (20, 300), + (22, 400), + ] + qft_configs = [(4,), (6,), (8,), (10,), (12,), (14,), (16,)] + + simulators = available_simulators(verbose=True) + + run_benchmarks("Random", random_configs, generate_random_gates, simulators) + run_benchmarks("QFT", qft_configs, generate_qft_gates, simulators) + + +if __name__ == "__main__": + main() diff --git a/benchmarks/plot_benchmarks.py b/benchmarks/plot_benchmarks.py new file mode 100644 index 00000000..eec6b6f0 --- /dev/null +++ b/benchmarks/plot_benchmarks.py @@ -0,0 +1,181 @@ +#!/usr/bin/env python3 +# Copyright 2025 qBraid +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Generate benchmark plots comparing PyQASM against other simulators. + +Third-party simulators (qiskit + qiskit-aer, cirq, pennylane) are optional; +any that are not installed are skipped and the plots show only the available +backends. Requires matplotlib. + +Usage: + python benchmarks/plot_benchmarks.py [--quick] [--out-dir DIR] +""" + +import argparse +import os +import sys + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +# Ensure we can import from the benchmarks directory +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) + +from bench_simulator import ( + available_simulators, + generate_qft_gates, + generate_random_gates, +) + +RANDOM_CONFIGS = [ + (4, 20), + (6, 40), + (8, 60), + (10, 80), + (12, 100), + (14, 150), + (16, 200), + (18, 250), + (20, 300), + (22, 400), +] + +QFT_CONFIGS = [4, 6, 8, 10, 12, 14, 16, 18, 20] + +QUICK_RANDOM_CONFIGS = [(2, 10), (4, 20), (6, 30)] +QUICK_QFT_CONFIGS = [2, 4, 6] + +COLORS = { + "PyQASM": "#2563eb", + "Qiskit Aer": "#16a34a", + "Cirq": "#eab308", + "PennyLane Lightning": "#dc2626", +} + +MARKERS = { + "PyQASM": "o", + "Qiskit Aer": "s", + "Cirq": "D", + "PennyLane Lightning": "^", +} + + +def collect_times(circuit_type, configs, simulators): + """Run benchmarks and return {simulator_name: ([qubits], [times_ms])}.""" + results = {name: ([], []) for name in simulators} + + for config in configs: + if circuit_type == "random": + nq, depth = config + num_qubits, gates = generate_random_gates(nq, depth) + label = f"n={nq}, d={depth}" + else: + nq = config + num_qubits, gates = generate_qft_gates(nq) + label = f"QFT n={nq}" + + print(f" {label} ...", end="", flush=True) + for name, (prepare_fn, bench_fn, _big_endian) in simulators.items(): + args = prepare_fn(num_qubits, gates) + t, _ = bench_fn(*args) + results[name][0].append(num_qubits) + results[name][1].append(t * 1000) # ms + print(f" {name}={t*1000:.1f}ms", end="", flush=True) + print() + + return results + + +def plot_comparison(results, title, out_path): + """Create a log-scale line plot comparing simulators.""" + fig, ax = plt.subplots(figsize=(10, 6)) + + for name, (qubits, times) in results.items(): + ax.plot( + qubits, + times, + marker=MARKERS.get(name, "o"), + color=COLORS.get(name), + linewidth=2, + markersize=7, + label=name, + ) + + ax.set_yscale("log") + ax.set_xlabel("Number of Qubits", fontsize=13) + ax.set_ylabel("Simulation Time (ms)", fontsize=13) + ax.set_title(title, fontsize=15) + ax.legend(fontsize=11) + ax.grid(True, which="both", ls="--", alpha=0.4) + ax.set_xticks(results["PyQASM"][0]) + fig.tight_layout() + fig.savefig(out_path, dpi=150) + print(f" Saved: {out_path}") + plt.close(fig) + + +def main(): + parser = argparse.ArgumentParser(description="Plot PyQASM simulator benchmarks") + parser.add_argument( + "--quick", + action="store_true", + help="Small qubit counts and fewer repeats (smoke test).", + ) + parser.add_argument( + "--out-dir", + default=os.path.dirname(os.path.abspath(__file__)), + help="Directory for the output PNGs (default: this benchmarks/ directory).", + ) + args = parser.parse_args() + + os.makedirs(args.out_dir, exist_ok=True) + + import bench_simulator + + if args.quick: + bench_simulator.N_REPEATS = 2 + random_configs = QUICK_RANDOM_CONFIGS + qft_configs = QUICK_QFT_CONFIGS + else: + random_configs = RANDOM_CONFIGS + qft_configs = QFT_CONFIGS + + simulators = available_simulators(verbose=True) + + print(f"Running benchmarks (median of {bench_simulator.N_REPEATS} runs each)\n") + + print("Random circuits:") + random_results = collect_times("random", random_configs, simulators) + plot_comparison( + random_results, + "Statevector Simulation Time — Random Circuits", + os.path.join(args.out_dir, "bench_random.png"), + ) + + print("\nQFT circuits:") + qft_results = collect_times("qft", qft_configs, simulators) + plot_comparison( + qft_results, + "Statevector Simulation Time — QFT Circuits", + os.path.join(args.out_dir, "bench_qft.png"), + ) + + print("\nDone.") + + +if __name__ == "__main__": + main() diff --git a/benchmarks/profile_simulator.py b/benchmarks/profile_simulator.py new file mode 100644 index 00000000..9497db0f --- /dev/null +++ b/benchmarks/profile_simulator.py @@ -0,0 +1,224 @@ +#!/usr/bin/env python3 +# Copyright 2025 qBraid +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Profiling script for the PyQASM statevector simulator. + +Three modes: + 1. cProfile function-level profiling of ``Simulator.run`` + 2. Section timing (preprocess vs simulation vs post-processing). Gate + fusion happens inline inside ``_preprocess``, so it is included in the + preprocess figure rather than reported separately. + 3. Per-kernel timing (accumulate time per opcode type) + +Usage: + python benchmarks/profile_simulator.py [--mode cprofile|section|kernel] \ + [--qubits N] [--depth D] +""" + +import argparse +import cProfile +import pstats +import random +import time +from io import StringIO + +import numpy as np + +from pyqasm import loads as pyqasm_loads +from pyqasm.accelerate.sv_sim import ( + apply_circuit, + apply_controlled_diagonal_gate, + apply_controlled_gate, + apply_diagonal_gate, + apply_single_qubit_gate, + apply_two_qubit_gate, +) +from pyqasm.simulator.statevector import ( + _OP_CONTROLLED, + _OP_CTRL_DIAGONAL, + _OP_DIAGONAL, + _OP_SINGLE, + _OP_TWO_QUBIT, + Simulator, + _preprocess, +) + +SINGLE_QUBIT_GATES = ["h", "x", "y", "z", "s", "t", "rx(1.0)", "ry(0.5)", "rz(0.3)"] +TWO_QUBIT_GATES = ["cx", "cy", "cz", "swap"] + +OP_NAMES = { + _OP_SINGLE: "single", + _OP_CONTROLLED: "controlled", + _OP_DIAGONAL: "diagonal", + _OP_CTRL_DIAGONAL: "ctrl_diag", + _OP_TWO_QUBIT: "two_qubit", +} + + +def generate_random_qasm(num_qubits: int, depth: int, seed: int = 42) -> str: + rng = random.Random(seed) + lines = [ + "OPENQASM 3;", + 'include "stdgates.inc";', + f"qubit[{num_qubits}] q;", + ] + for _ in range(depth): + if num_qubits >= 2 and rng.random() < 0.4: + gate = rng.choice(TWO_QUBIT_GATES) + q0, q1 = rng.sample(range(num_qubits), 2) + lines.append(f"{gate} q[{q0}], q[{q1}];") + else: + gate = rng.choice(SINGLE_QUBIT_GATES) + q = rng.randint(0, num_qubits - 1) + lines.append(f"{gate} q[{q}];") + return "\n".join(lines) + + +def _prepared_module(num_qubits: int, depth: int): + """Load, unroll, and strip idle qubits from a random circuit.""" + qasm = generate_random_qasm(num_qubits, depth) + module = pyqasm_loads(qasm) + module.unroll() + module.remove_idle_qubits() + return module + + +def mode_cprofile(num_qubits, depth): + """Function-level profiling with cProfile.""" + module = _prepared_module(num_qubits, depth) + sim = Simulator(seed=0) + + print(f"\n=== cProfile: {num_qubits} qubits, depth {depth} ===\n") + pr = cProfile.Profile() + pr.enable() + sim.run(module, shots=0) + pr.disable() + + s = StringIO() + ps = pstats.Stats(pr, stream=s).sort_stats("tottime") + ps.print_stats(30) + print(s.getvalue()) + + +def mode_section(num_qubits, depth): + """Section timing: preprocess (incl. inline fusion) vs simulation vs post.""" + module = _prepared_module(num_qubits, depth) + nq = module.num_qubits + + print(f"\n=== Section timing: {num_qubits} qubits, depth {depth} ===\n") + + # Preprocess (gate fusion happens inline here) + t0 = time.perf_counter_ns() + n, opcodes, targets, controls, gate_params, diag_phases, tq_offsets, tq_gates = _preprocess( + module, nq + ) + t1 = time.perf_counter_ns() + + # Simulation + sv = np.zeros(2**nq, dtype=np.complex128) + sv[0] = 1.0 + if n > 0: + apply_circuit( + sv, nq, opcodes, targets, controls, gate_params, diag_phases, tq_offsets, tq_gates, n + ) + t2 = time.perf_counter_ns() + + # Post-processing + probabilities = np.abs(sv) ** 2 # noqa: F841 (timed for its cost, not its value) + t3 = time.perf_counter_ns() + + pre_us = (t1 - t0) / 1000 + sim_us = (t2 - t1) / 1000 + post_us = (t3 - t2) / 1000 + total_us = (t3 - t0) / 1000 + + print(f" Instructions: {n}") + print(f" Preprocess: {pre_us:10.1f} us ({pre_us/total_us*100:5.1f}%)") + print(f" Simulation: {sim_us:10.1f} us ({sim_us/total_us*100:5.1f}%)") + print(f" Post-process: {post_us:10.1f} us ({post_us/total_us*100:5.1f}%)") + print(f" Total: {total_us:10.1f} us") + + +def mode_kernel(num_qubits, depth): + """Per-kernel timing: accumulate time per opcode type.""" + module = _prepared_module(num_qubits, depth) + nq = module.num_qubits + + n, opcodes, targets, controls, gate_params, diag_phases, tq_offsets, tq_gates = _preprocess( + module, nq + ) + + sv = np.zeros(2**nq, dtype=np.complex128) + sv[0] = 1.0 + + kernel_times = {} + kernel_counts = {} + + print(f"\n=== Per-kernel timing: {num_qubits} qubits, depth {depth} ===\n") + + for i in range(n): + op = int(opcodes[i]) + tgt = int(targets[i]) + ctrl = int(controls[i]) + gp_off = i * 4 + dp_off = i * 2 + + t0 = time.perf_counter_ns() + if op == _OP_SINGLE: + flat = gate_params[gp_off : gp_off + 4].copy() + apply_single_qubit_gate(sv, nq, tgt, flat) + elif op == _OP_CONTROLLED: + flat = gate_params[gp_off : gp_off + 4].copy() + apply_controlled_gate(sv, nq, ctrl, tgt, flat) + elif op == _OP_DIAGONAL: + apply_diagonal_gate(sv, nq, tgt, diag_phases[dp_off], diag_phases[dp_off + 1]) + elif op == _OP_CTRL_DIAGONAL: + apply_controlled_diagonal_gate(sv, nq, ctrl, tgt, diag_phases[dp_off]) + elif op == _OP_TWO_QUBIT: + tq_off = int(tq_offsets[i]) + flat = tq_gates[tq_off : tq_off + 16].copy() + apply_two_qubit_gate(sv, nq, ctrl, tgt, flat) + t1 = time.perf_counter_ns() + + name = OP_NAMES.get(op, f"op_{op}") + kernel_times[name] = kernel_times.get(name, 0) + (t1 - t0) + kernel_counts[name] = kernel_counts.get(name, 0) + 1 + + total_ns = sum(kernel_times.values()) + for name in sorted(kernel_times.keys()): + t_us = kernel_times[name] / 1000 + cnt = kernel_counts[name] + pct = kernel_times[name] / total_ns * 100 if total_ns > 0 else 0 + print(f" {name:15s}: {t_us:10.1f} us ({pct:5.1f}%) count={cnt}") + print(f" {'total':15s}: {total_ns/1000:10.1f} us") + + +def main(): + parser = argparse.ArgumentParser(description="Profile PyQASM simulator") + parser.add_argument("--mode", choices=["cprofile", "section", "kernel"], default="section") + parser.add_argument("--qubits", type=int, default=16) + parser.add_argument("--depth", type=int, default=200) + args = parser.parse_args() + + if args.mode == "cprofile": + mode_cprofile(args.qubits, args.depth) + elif args.mode == "section": + mode_section(args.qubits, args.depth) + elif args.mode == "kernel": + mode_kernel(args.qubits, args.depth) + + +if __name__ == "__main__": + main() diff --git a/bin/cibw/build_wheels.sh b/bin/cibw/build_wheels.sh index 5accab20..1b61e97a 100755 --- a/bin/cibw/build_wheels.sh +++ b/bin/cibw/build_wheels.sh @@ -18,7 +18,9 @@ set -e echo "Running build_wheels.sh" -python -m pip install cibuildwheel +# Pin the major version: cibuildwheel sits in the release critical path and +# its defaults (images, repair commands) change across majors. +python -m pip install "cibuildwheel>=3.0,<4" python -m cibuildwheel --output-dir dist build_exit_code=$? diff --git a/bin/cibw/test_wheel.sh b/bin/cibw/test_wheel.sh index e1fb091f..f649398f 100755 --- a/bin/cibw/test_wheel.sh +++ b/bin/cibw/test_wheel.sh @@ -28,8 +28,12 @@ echo "Running test_wheel.sh" project=$1 # the built wheel is already installed in the test env by cibuildwheel -# just run the tests -python -m pytest $project/tests +# just run the tests. +# Benchmark/perf-regression tests enforce hardware-dependent time budgets and are +# excluded from CI (see tests/test_perf_regression.py). The qiskit comparison tests +# skip automatically when qiskit is unavailable (it is not installable in the +# manylinux wheel-build containers). +python -m pytest "$project/tests" -m "not benchmark" pytest_exit_code=$? if [ $pytest_exit_code -ne 0 ]; then diff --git a/bin/test_sdist.sh b/bin/test_sdist.sh index eaaa5ae6..5824794d 100755 --- a/bin/test_sdist.sh +++ b/bin/test_sdist.sh @@ -32,7 +32,7 @@ source "$TEMP_ENV_DIR/bin/activate" # Install the source distribution SCRIPT_DIR="$TARGET_PATH/bin" -"$SCRIPT_DIR/install_wheel_extras.sh" "$TARGET_PATH/dist" --type sdist --extra cli --extra test --extra pulse +"$SCRIPT_DIR/install_wheel_extras.sh" "$TARGET_PATH/dist" --type sdist --extra cli --extra test --extra test-sim --extra pulse # Print the installed version python -c "import pyqasm; print('Installed pyqasm version:', pyqasm.__version__)" @@ -56,8 +56,10 @@ if [[ ${RELEASE_BUILD:-false} == "true" ]]; then fi fi -# Run the tests on the installed source distribution -pytest "$TARGET_PATH/tests" +# Run the tests on the installed source distribution. +# Benchmark/perf-regression tests enforce hardware-dependent time budgets and are +# excluded from CI (see tests/test_perf_regression.py). +pytest "$TARGET_PATH/tests" -m "not benchmark" # Deactivate and remove venv deactivate diff --git a/docs/api/pyqasm.simulator.rst b/docs/api/pyqasm.simulator.rst new file mode 100644 index 00000000..e0ea5955 --- /dev/null +++ b/docs/api/pyqasm.simulator.rst @@ -0,0 +1,4 @@ +pyqasm.simulator +================= + +.. automodule:: pyqasm.simulator diff --git a/docs/index.rst b/docs/index.rst index 766bffb7..320970bf 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -192,6 +192,7 @@ Resources :hidden: api/pyqasm + api/pyqasm.simulator .. toctree:: :caption: ALGOS API Reference diff --git a/pyproject.toml b/pyproject.toml index 88926526..5f10accf 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -43,19 +43,22 @@ Homepage = "https://docs.qbraid.com/v2/pyqasm" [project.optional-dependencies] cli = ["typer>=0.12.1", "rich>=10.11.0", "typing-extensions"] test = ["pytest", "pytest-cov", "pytest-mpl", "pillow<12.3.0", "matplotlib", "tabulate"] +# Reference-oracle deps used only by the statevector simulator comparison tests. +# Kept separate from `test` because qiskit/qiskit-aer have no prebuilt wheels for +# the manylinux wheel-build containers and would require a Rust toolchain there. +test-sim = ["qiskit", "qiskit-aer", "qiskit-qasm3-import"] lint = ["black", "isort>=6.0.0", "pylint", "mypy", "qbraid-cli>=0.10.2"] docs = ["sphinx>=7.3.7,<8.3.0", "sphinx-autodoc-typehints>=1.24,<3.2", "sphinx-rtd-theme>=2.0.0,<4.0.0", "docutils<0.23", "sphinx-copybutton"] +simulation = ["numba>=0.59"] visualization = ["pillow<12.3.0", "matplotlib", "tabulate"] pulse = ["openpulse>=1.0.1"] [tool.setuptools.package-data] pyqasm = ["py.typed", "*.pyx"] +"pyqasm.accelerate" = ["*.pyx"] [tool.setuptools] packages = {find = {where = ["src"]}} -ext-modules = [ - {name = "pyqasm.accelerate.linalg", sources = ["src/pyqasm/accelerate/linalg.pyx"], include-dirs = ["numpy"]} -] [project.scripts] pyqasm = "pyqasm.cli.main:app" @@ -75,6 +78,9 @@ line_length = 100 [tool.pytest.ini_options] addopts = "-ra" testpaths = ["tests"] +markers = [ + "benchmark: performance regression tests (may be slow)", +] [tool.coverage.run] source = ["pyqasm"] diff --git a/setup.py b/setup.py new file mode 100644 index 00000000..d8716df4 --- /dev/null +++ b/setup.py @@ -0,0 +1,181 @@ +# Copyright 2025 qBraid +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Build configuration for Cython extensions with optional OpenMP support. + +Optimization flags are chosen for portability of the distributed wheels: + + * No ``-march`` (and certainly not ``-march=native``). Wheels are built at the + baseline ISA so they run on any CPU of the target architecture, matching what + SciPy and scikit-learn do. ``-march=native`` would tie a wheel to the exact + CPU of the build machine and can ``SIGILL`` on older hardware. + * No ``-ffast-math``. The statevector simulator relies on IEEE-754 semantics + (NaN/Inf handling, no float reassociation), so we keep strict FP. + +Flags are also selected per-compiler: ``-O3``/``-fopenmp`` are GCC/Clang +spellings and are silently ignored by MSVC, which needs ``/O2``/``/openmp``. + +OpenMP policy per platform: + + * Linux: enabled when the probe succeeds (the norm on manylinux). + * macOS: disabled by default. There is no portable system libomp, and + bundling Homebrew's breaks ``delocate``. Set ``PYQASM_MACOS_OPENMP=1`` for + a local, non-distributed build. + * Windows: disabled by default. MSVC's ``/openmp`` makes the extension + depend on ``vcomp140.dll``, which is not part of a stock Windows install; + ``delvewheel`` then vendors whichever copy it finds first on PATH (on + GitHub runners that has been ImageMagick's), making releases + nondeterministic. Set ``PYQASM_WINDOWS_OPENMP=1`` to opt in locally. + +OpenMP only affects compiler flags: the generated C guards all OpenMP pragmas +behind ``#ifdef _OPENMP``, so the same sources build serial or parallel and no +flavor is baked into sdists. +""" + +import os +import shlex +import subprocess +import sys +import tempfile + +import numpy as np +from setuptools import Extension, setup +from setuptools.command.build_ext import build_ext + +# Extensions whose kernels are OpenMP-parallelized. +_OPENMP_EXTENSIONS = {"pyqasm.accelerate.sv_sim"} + + +def _detect_openmp_unix() -> tuple[list[str], list[str]]: + """Detect OpenMP availability for GCC/Clang. + + Takes no parameters; probes the compiler named by ``$CC`` (or ``cc``) with + candidate flag sets. Returns ``(compile_args, link_args)`` for the first + candidate that compiles a test program, or ``([], [])`` when OpenMP is + unavailable or disabled by platform policy. + """ + if sys.platform == "darwin" and not os.environ.get("PYQASM_MACOS_OPENMP"): + return [], [] + + test_code = b"#include \nint main() { return omp_get_max_threads(); }\n" + + # Probe with the compiler the build will actually use: honor $CC (which may + # include flags, e.g. "gcc -pthread") before falling back to plain `cc`. + compiler_cmd = shlex.split(os.environ.get("CC") or "cc") + + candidates = [] + if sys.platform == "darwin": + # Opt-in macOS OpenMP: libomp from Homebrew. + for prefix in ("/opt/homebrew/opt/libomp", "/usr/local/opt/libomp"): + if os.path.isdir(prefix): + candidates.append( + ( + ["-Xpreprocessor", "-fopenmp", f"-I{prefix}/include"], + [f"-L{prefix}/lib", "-lomp"], + ) + ) + # Also try Xcode / system clang with -fopenmp + candidates.append((["-fopenmp"], ["-fopenmp"])) + else: + # Linux / other: standard -fopenmp + candidates.append((["-fopenmp"], ["-fopenmp"])) + + for cflags, ldflags in candidates: + tmp_name = None + try: + with tempfile.NamedTemporaryFile(suffix=".c", delete=False) as f: + f.write(test_code) + f.flush() + tmp_name = f.name + cmd = compiler_cmd + cflags + ldflags + [tmp_name, "-o", tmp_name + ".out"] + subprocess.check_call(cmd, stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL) + os.unlink(tmp_name + ".out") + return cflags, ldflags + except (subprocess.CalledProcessError, FileNotFoundError, OSError): + continue + finally: + if tmp_name and os.path.exists(tmp_name): + os.unlink(tmp_name) + + return [], [] + + +if sys.platform.startswith("win"): + # MSVC: /openmp is always available but opt-in (see module docstring). + USE_OPENMP = bool(os.environ.get("PYQASM_WINDOWS_OPENMP")) + _OMP_COMPILE_UNIX, _OMP_LINK_UNIX = [], [] +else: + _OMP_COMPILE_UNIX, _OMP_LINK_UNIX = _detect_openmp_unix() + USE_OPENMP = bool(_OMP_COMPILE_UNIX) + + +class BuildExt(build_ext): + """Cythonize lazily and inject portable, compiler-appropriate flags. + + Running ``cythonize`` here rather than at module scope keeps the ``.pyx`` + files in ``Extension.sources`` while sdists and wheels are being assembled, + so generated ``.c`` never ships in release artifacts, and metadata-only + PEP 517 hooks skip Cython codegen entirely. + """ + + def finalize_options(self) -> None: + """Cythonize the .pyx extensions in place, then finalize as usual.""" + # pylint: disable-next=import-outside-toplevel + from Cython.Build import cythonize + + self.distribution.ext_modules[:] = cythonize( + self.distribution.ext_modules, + language_level=3, + ) + super().finalize_options() + + def build_extensions(self) -> None: + """Apply per-compiler optimization and OpenMP flags, then build.""" + is_msvc = self.compiler.compiler_type == "msvc" + base_compile = ["/O2"] if is_msvc else ["-O3"] + if USE_OPENMP: + omp_compile, omp_link = ( + (["/openmp"], []) if is_msvc else (_OMP_COMPILE_UNIX, _OMP_LINK_UNIX) + ) + print(f"OpenMP enabled: compile={omp_compile}, link={omp_link}") + else: + omp_compile, omp_link = [], [] + print("OpenMP disabled - building single-threaded kernels") + + for ext in self.extensions: + ext.extra_compile_args = base_compile + list(ext.extra_compile_args) + if ext.name in _OPENMP_EXTENSIONS and USE_OPENMP: + ext.extra_compile_args += omp_compile + ext.extra_link_args = list(ext.extra_link_args) + omp_link + + super().build_extensions() + + +extensions = [ + Extension( + "pyqasm.accelerate.linalg", + sources=["src/pyqasm/accelerate/linalg.pyx"], + include_dirs=[np.get_include()], + ), + Extension( + "pyqasm.accelerate.sv_sim", + sources=["src/pyqasm/accelerate/sv_sim.pyx"], + include_dirs=[np.get_include()], + ), +] + +setup( + ext_modules=extensions, + cmdclass={"build_ext": BuildExt}, +) diff --git a/src/pyqasm/accelerate/sv_sim.pyx b/src/pyqasm/accelerate/sv_sim.pyx new file mode 100644 index 00000000..f4fb5175 --- /dev/null +++ b/src/pyqasm/accelerate/sv_sim.pyx @@ -0,0 +1,515 @@ +# Copyright 2025 qBraid +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# cython: language_level=3 +# cython: infer_types=True + +import os + +import cython + +from cython.parallel cimport prange + +# OpenMP support is decided purely by compiler flags (-fopenmp / /openmp) in +# setup.py. Cython emits the OpenMP pragmas for prange behind `#ifdef _OPENMP` +# guards, so the same generated C compiles either way: with OpenMP enabled the +# prange paths parallelize, without it they degrade to serial loops. Nothing +# here cimports the `openmp` module, so no unguarded is ever pulled in +# and no process-global OpenMP state is mutated at import time. + +# Parallelism is opt-in via the PYQASM_NUM_THREADS env var, read once at +# import. The default is 1 (serial) to avoid conflicts with other OpenMP +# runtimes (e.g. Qiskit Aer) in the same process. The count is stored in a +# module-level C int and passed explicitly to prange(num_threads=...), so the +# policy does not depend on any per-thread OpenMP internal control variable +# and behaves identically on every Python thread. + + +def _read_num_threads(): + """Parse PYQASM_NUM_THREADS defensively; malformed or non-positive -> 1.""" + raw = os.environ.get("PYQASM_NUM_THREADS", "1") + try: + value = int(raw) + except (TypeError, ValueError): + return 1 + return value if value > 0 else 1 + + +cdef int _NUM_THREADS = _read_num_threads() + +# Threshold: parallelize when the statevector has >= this many pairs. +# 2^17 = 131072 pairs corresponds to ~18 qubits (4 MB statevector). +# Below this, OpenMP thread overhead dominates over parallel gains. +cdef Py_ssize_t PARALLEL_THRESHOLD = 131072 + + +# --- cdef kernel functions --- +# Each kernel manages its own GIL release: prange paths use prange's +# built-in nogil, serial paths use explicit `with nogil:` blocks. +# Kernels are noexcept raw-pointer code; ALL argument validation happens in +# the cpdef wrappers below before any pointer is taken. + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef void _apply_single_qubit_gate( + double complex* sv, + Py_ssize_t num_qubits, + Py_ssize_t target, + double complex g00, + double complex g01, + double complex g10, + double complex g11, +) noexcept: + cdef Py_ssize_t half = 1 << (num_qubits - 1) + cdef Py_ssize_t step = 1 << target + cdef Py_ssize_t s, i0, i1 + cdef double complex a0, a1 + + if _NUM_THREADS > 1 and half >= PARALLEL_THRESHOLD: + for s in prange(half, schedule='static', num_threads=_NUM_THREADS, nogil=True): + i0 = ((s >> target) << (target + 1)) | (s & (step - 1)) + i1 = i0 | step + a0 = sv[i0] + a1 = sv[i1] + sv[i0] = g00 * a0 + g01 * a1 + sv[i1] = g10 * a0 + g11 * a1 + return + with nogil: + for s in range(half): + i0 = ((s >> target) << (target + 1)) | (s & (step - 1)) + i1 = i0 | step + a0 = sv[i0] + a1 = sv[i1] + sv[i0] = g00 * a0 + g01 * a1 + sv[i1] = g10 * a0 + g11 * a1 + + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef void _apply_controlled_gate( + double complex* sv, + Py_ssize_t num_qubits, + Py_ssize_t control, + Py_ssize_t target, + double complex g00, + double complex g01, + double complex g10, + double complex g11, +) noexcept: + cdef Py_ssize_t quarter = 1 << (num_qubits - 2) + cdef Py_ssize_t lo, hi + cdef Py_ssize_t lo_mask, hi_limit + cdef Py_ssize_t s, base, i0, i1 + cdef double complex a0, a1 + cdef Py_ssize_t ctrl_bit = 1 << control + cdef Py_ssize_t tgt_bit = 1 << target + + if control < target: + lo = control + hi = target + else: + lo = target + hi = control + + lo_mask = (1 << lo) - 1 + hi_limit = (1 << hi) - 1 + + if _NUM_THREADS > 1 and quarter >= PARALLEL_THRESHOLD: + for s in prange(quarter, schedule='static', num_threads=_NUM_THREADS, nogil=True): + base = (s & lo_mask) | (((s >> lo) << (lo + 1)) & hi_limit) | ((s >> (hi - 1)) << (hi + 1)) + base = base | ctrl_bit + i0 = base + i1 = base | tgt_bit + a0 = sv[i0] + a1 = sv[i1] + sv[i0] = g00 * a0 + g01 * a1 + sv[i1] = g10 * a0 + g11 * a1 + return + with nogil: + for s in range(quarter): + base = (s & lo_mask) | (((s >> lo) << (lo + 1)) & hi_limit) | ((s >> (hi - 1)) << (hi + 1)) + base = base | ctrl_bit + i0 = base + i1 = base | tgt_bit + a0 = sv[i0] + a1 = sv[i1] + sv[i0] = g00 * a0 + g01 * a1 + sv[i1] = g10 * a0 + g11 * a1 + + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef void _apply_diagonal_gate( + double complex* sv, + Py_ssize_t num_qubits, + Py_ssize_t target, + double complex phase0, + double complex phase1, +) noexcept: + cdef Py_ssize_t half = 1 << (num_qubits - 1) + cdef Py_ssize_t step = 1 << target + cdef Py_ssize_t s, i0, i1 + + if _NUM_THREADS > 1 and half >= PARALLEL_THRESHOLD: + for s in prange(half, schedule='static', num_threads=_NUM_THREADS, nogil=True): + i0 = ((s >> target) << (target + 1)) | (s & (step - 1)) + i1 = i0 | step + sv[i0] = phase0 * sv[i0] + sv[i1] = phase1 * sv[i1] + return + with nogil: + for s in range(half): + i0 = ((s >> target) << (target + 1)) | (s & (step - 1)) + i1 = i0 | step + sv[i0] = phase0 * sv[i0] + sv[i1] = phase1 * sv[i1] + + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef void _apply_controlled_diagonal_gate( + double complex* sv, + Py_ssize_t num_qubits, + Py_ssize_t control, + Py_ssize_t target, + double complex phase, +) noexcept: + cdef Py_ssize_t quarter = 1 << (num_qubits - 2) + cdef Py_ssize_t lo, hi + cdef Py_ssize_t lo_mask, hi_limit + cdef Py_ssize_t s, base, idx + cdef Py_ssize_t ctrl_bit = 1 << control + cdef Py_ssize_t tgt_bit = 1 << target + + if control < target: + lo = control + hi = target + else: + lo = target + hi = control + + lo_mask = (1 << lo) - 1 + hi_limit = (1 << hi) - 1 + + if _NUM_THREADS > 1 and quarter >= PARALLEL_THRESHOLD: + for s in prange(quarter, schedule='static', num_threads=_NUM_THREADS, nogil=True): + base = (s & lo_mask) | (((s >> lo) << (lo + 1)) & hi_limit) | ((s >> (hi - 1)) << (hi + 1)) + idx = base | ctrl_bit | tgt_bit + sv[idx] = phase * sv[idx] + return + with nogil: + for s in range(quarter): + base = (s & lo_mask) | (((s >> lo) << (lo + 1)) & hi_limit) | ((s >> (hi - 1)) << (hi + 1)) + idx = base | ctrl_bit | tgt_bit + sv[idx] = phase * sv[idx] + + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef void _apply_two_qubit_gate( + double complex* sv, + Py_ssize_t num_qubits, + Py_ssize_t qubit0, + Py_ssize_t qubit1, + double complex* gate, +) noexcept: + cdef Py_ssize_t quarter = 1 << (num_qubits - 2) + cdef Py_ssize_t lo, hi + cdef Py_ssize_t lo_mask, hi_limit + cdef Py_ssize_t s, base, i00, i01, i10, i11 + cdef double complex a00, a01, a10, a11 + cdef Py_ssize_t q0_bit = 1 << qubit0 + cdef Py_ssize_t q1_bit = 1 << qubit1 + cdef double complex g0 = gate[0], g1 = gate[1], g2 = gate[2], g3 = gate[3] + cdef double complex g4 = gate[4], g5 = gate[5], g6 = gate[6], g7 = gate[7] + cdef double complex g8 = gate[8], g9 = gate[9], g10 = gate[10], g11_ = gate[11] + cdef double complex g12 = gate[12], g13 = gate[13], g14 = gate[14], g15 = gate[15] + + if qubit0 < qubit1: + lo = qubit0 + hi = qubit1 + else: + lo = qubit1 + hi = qubit0 + + lo_mask = (1 << lo) - 1 + hi_limit = (1 << hi) - 1 + + if _NUM_THREADS > 1 and quarter >= PARALLEL_THRESHOLD: + for s in prange(quarter, schedule='static', num_threads=_NUM_THREADS, nogil=True): + base = (s & lo_mask) | (((s >> lo) << (lo + 1)) & hi_limit) | ((s >> (hi - 1)) << (hi + 1)) + i00 = base + i01 = base | q0_bit + i10 = base | q1_bit + i11 = base | q0_bit | q1_bit + a00 = sv[i00] + a01 = sv[i01] + a10 = sv[i10] + a11 = sv[i11] + sv[i00] = g0 * a00 + g1 * a01 + g2 * a10 + g3 * a11 + sv[i01] = g4 * a00 + g5 * a01 + g6 * a10 + g7 * a11 + sv[i10] = g8 * a00 + g9 * a01 + g10 * a10 + g11_ * a11 + sv[i11] = g12 * a00 + g13 * a01 + g14 * a10 + g15 * a11 + return + with nogil: + for s in range(quarter): + base = (s & lo_mask) | (((s >> lo) << (lo + 1)) & hi_limit) | ((s >> (hi - 1)) << (hi + 1)) + i00 = base + i01 = base | q0_bit + i10 = base | q1_bit + i11 = base | q0_bit | q1_bit + a00 = sv[i00] + a01 = sv[i01] + a10 = sv[i10] + a11 = sv[i11] + sv[i00] = g0 * a00 + g1 * a01 + g2 * a10 + g3 * a11 + sv[i01] = g4 * a00 + g5 * a01 + g6 * a10 + g7 * a11 + sv[i10] = g8 * a00 + g9 * a01 + g10 * a10 + g11_ * a11 + sv[i11] = g12 * a00 + g13 * a01 + g14 * a10 + g15 * a11 + + +# --- argument validation --- +# The kernels above are unchecked raw-pointer code, so every cpdef entry point +# validates its arguments first: without this, a malformed call from Python +# segfaults the interpreter instead of raising. + +cdef inline void _check_state(double complex[::1] sv, Py_ssize_t num_qubits) except *: + if num_qubits < 1 or num_qubits > 62: + raise ValueError(f"num_qubits must be in [1, 62], got {num_qubits}.") + if sv.shape[0] != ( 1) << num_qubits: + raise ValueError( + f"Statevector length {sv.shape[0]} does not match " + f"2**num_qubits = {( 1) << num_qubits}." + ) + + +cdef inline void _check_target(Py_ssize_t num_qubits, Py_ssize_t target) except *: + if target < 0 or target >= num_qubits: + raise ValueError( + f"Target qubit {target} out of range for {num_qubits} qubit(s)." + ) + + +cdef inline void _check_control_target( + Py_ssize_t num_qubits, Py_ssize_t control, Py_ssize_t target +) except *: + _check_target(num_qubits, target) + if control < 0 or control >= num_qubits: + raise ValueError( + f"Control qubit {control} out of range for {num_qubits} qubit(s)." + ) + if control == target: + raise ValueError(f"Control and target qubit are both {target}.") + + +# --- cpdef wrappers for Python API --- + +cpdef void apply_single_qubit_gate( + double complex[::1] sv, + Py_ssize_t num_qubits, + Py_ssize_t target, + double complex[::1] gate, +): + _check_state(sv, num_qubits) + _check_target(num_qubits, target) + if gate.shape[0] < 4: + raise ValueError(f"Gate buffer must hold 4 elements, got {gate.shape[0]}.") + cdef double complex* sv_ptr = &sv[0] + _apply_single_qubit_gate(sv_ptr, num_qubits, target, gate[0], gate[1], gate[2], gate[3]) + + +cpdef void apply_controlled_gate( + double complex[::1] sv, + Py_ssize_t num_qubits, + Py_ssize_t control, + Py_ssize_t target, + double complex[::1] gate, +): + _check_state(sv, num_qubits) + _check_control_target(num_qubits, control, target) + if gate.shape[0] < 4: + raise ValueError(f"Gate buffer must hold 4 elements, got {gate.shape[0]}.") + cdef double complex* sv_ptr = &sv[0] + _apply_controlled_gate(sv_ptr, num_qubits, control, target, gate[0], gate[1], gate[2], gate[3]) + + +cpdef void apply_diagonal_gate( + double complex[::1] sv, + Py_ssize_t num_qubits, + Py_ssize_t target, + double complex phase0, + double complex phase1, +): + _check_state(sv, num_qubits) + _check_target(num_qubits, target) + cdef double complex* sv_ptr = &sv[0] + _apply_diagonal_gate(sv_ptr, num_qubits, target, phase0, phase1) + + +cpdef void apply_controlled_diagonal_gate( + double complex[::1] sv, + Py_ssize_t num_qubits, + Py_ssize_t control, + Py_ssize_t target, + double complex phase, +): + _check_state(sv, num_qubits) + _check_control_target(num_qubits, control, target) + cdef double complex* sv_ptr = &sv[0] + _apply_controlled_diagonal_gate(sv_ptr, num_qubits, control, target, phase) + + +cpdef void apply_two_qubit_gate( + double complex[::1] sv, + Py_ssize_t num_qubits, + Py_ssize_t qubit0, + Py_ssize_t qubit1, + double complex[::1] gate, +): + _check_state(sv, num_qubits) + _check_control_target(num_qubits, qubit0, qubit1) + if gate.shape[0] < 16: + raise ValueError(f"Gate buffer must hold 16 elements, got {gate.shape[0]}.") + cdef double complex* sv_ptr = &sv[0] + cdef double complex* gate_ptr = &gate[0] + _apply_two_qubit_gate(sv_ptr, num_qubits, qubit0, qubit1, gate_ptr) + + +# --- Batch dispatch: single Python->C crossing for entire circuit --- + +cdef enum OpCode: + OP_SINGLE = 0 + OP_CONTROLLED = 1 + OP_DIAGONAL = 2 + OP_CTRL_DIAGONAL = 3 + OP_TWO_QUBIT = 4 + + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cpdef void apply_circuit( + double complex[::1] sv, + Py_ssize_t num_qubits, + int[::1] opcodes, + int[::1] targets, + int[::1] controls, + double complex[::1] gate_params, + double complex[::1] diag_phases, + int[::1] two_qubit_offsets, + double complex[::1] two_qubit_gates, + Py_ssize_t n_instructions, +): + """Execute an entire circuit. Each gate releases the GIL internally. + + All instructions are validated up front (bounds, opcode range, buffer + sizes) before any kernel touches the statevector, so invalid input raises + ValueError instead of corrupting memory. + """ + cdef Py_ssize_t i, gp_offset, dp_offset, tq_offset + cdef int op, tgt, ctrl + + # --- validation pass (GIL held, O(n_instructions)) --- + _check_state(sv, num_qubits) + if n_instructions < 0: + raise ValueError(f"n_instructions must be >= 0, got {n_instructions}.") + if ( + opcodes.shape[0] < n_instructions + or targets.shape[0] < n_instructions + or controls.shape[0] < n_instructions + or two_qubit_offsets.shape[0] < n_instructions + ): + raise ValueError( + f"Instruction arrays are shorter than n_instructions = {n_instructions}." + ) + if gate_params.shape[0] < n_instructions * 4 or diag_phases.shape[0] < n_instructions * 2: + raise ValueError( + f"Parameter arrays are shorter than required for " + f"n_instructions = {n_instructions}." + ) + for i in range(n_instructions): + op = opcodes[i] + tgt = targets[i] + ctrl = controls[i] + if op < OP_SINGLE or op > OP_TWO_QUBIT: + raise ValueError(f"Instruction {i}: invalid opcode {op}.") + if tgt < 0 or tgt >= num_qubits: + raise ValueError( + f"Instruction {i}: target qubit {tgt} out of range for " + f"{num_qubits} qubit(s)." + ) + if op == OP_CONTROLLED or op == OP_CTRL_DIAGONAL or op == OP_TWO_QUBIT: + if ctrl < 0 or ctrl >= num_qubits: + raise ValueError( + f"Instruction {i}: control qubit {ctrl} out of range for " + f"{num_qubits} qubit(s)." + ) + if ctrl == tgt: + raise ValueError( + f"Instruction {i}: control and target qubit are both {tgt}." + ) + if op == OP_TWO_QUBIT: + tq_offset = two_qubit_offsets[i] + if tq_offset < 0 or tq_offset + 16 > two_qubit_gates.shape[0]: + raise ValueError( + f"Instruction {i}: two-qubit gate offset {tq_offset} out of " + f"range for buffer of length {two_qubit_gates.shape[0]}." + ) + + # --- dispatch pass --- + cdef double complex* sv_ptr = &sv[0] + cdef double complex* gp_ptr = &gate_params[0] + cdef double complex* dp_ptr = &diag_phases[0] + cdef double complex* tq_ptr = &two_qubit_gates[0] + + for i in range(n_instructions): + op = opcodes[i] + tgt = targets[i] + ctrl = controls[i] + gp_offset = i * 4 + dp_offset = i * 2 + + if op == OP_SINGLE: + _apply_single_qubit_gate( + sv_ptr, num_qubits, tgt, + gp_ptr[gp_offset], gp_ptr[gp_offset + 1], + gp_ptr[gp_offset + 2], gp_ptr[gp_offset + 3], + ) + elif op == OP_CONTROLLED: + _apply_controlled_gate( + sv_ptr, num_qubits, ctrl, tgt, + gp_ptr[gp_offset], gp_ptr[gp_offset + 1], + gp_ptr[gp_offset + 2], gp_ptr[gp_offset + 3], + ) + elif op == OP_DIAGONAL: + _apply_diagonal_gate( + sv_ptr, num_qubits, tgt, + dp_ptr[dp_offset], dp_ptr[dp_offset + 1], + ) + elif op == OP_CTRL_DIAGONAL: + _apply_controlled_diagonal_gate( + sv_ptr, num_qubits, ctrl, tgt, + dp_ptr[dp_offset], + ) + elif op == OP_TWO_QUBIT: + tq_offset = two_qubit_offsets[i] + _apply_two_qubit_gate( + sv_ptr, num_qubits, ctrl, tgt, + &tq_ptr[tq_offset], + ) diff --git a/src/pyqasm/simulator/__init__.py b/src/pyqasm/simulator/__init__.py new file mode 100644 index 00000000..a1d56354 --- /dev/null +++ b/src/pyqasm/simulator/__init__.py @@ -0,0 +1,33 @@ +# Copyright 2026 qBraid +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +""" +Module containing the PyQASM statevector simulator. + +.. currentmodule:: pyqasm.simulator + +Classes +--------- + +.. autosummary:: + :toctree: ../stubs/ + + Simulator + SimulatorResult + +""" + +from .statevector import Simulator, SimulatorResult + +__all__ = ["Simulator", "SimulatorResult"] diff --git a/src/pyqasm/simulator/statevector.py b/src/pyqasm/simulator/statevector.py new file mode 100644 index 00000000..eb355aff --- /dev/null +++ b/src/pyqasm/simulator/statevector.py @@ -0,0 +1,818 @@ +# Copyright 2025 qBraid +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +""" +Statevector simulator for PyQASM. + +""" + +# pylint: disable=no-name-in-module,too-many-return-statements +# pylint: disable=too-many-locals,too-many-branches,too-many-statements +# pylint: disable=arguments-out-of-order,unused-argument + +import os +from collections import Counter +from dataclasses import dataclass + +import numpy as np +from openqasm3.ast import ( + BinaryExpression, + BinaryOperator, + BranchingStatement, + ClassicalDeclaration, + FloatLiteral, + Identifier, + Include, + IndexedIdentifier, + IntegerLiteral, + QuantumBarrier, + QuantumGate, + QuantumMeasurementStatement, + QuantumReset, + QubitDeclaration, + UnaryExpression, + UnaryOperator, +) + +from pyqasm import loads +from pyqasm.accelerate.sv_sim import apply_circuit # type: ignore[import-not-found] +from pyqasm.modules.base import QasmModule + +try: + import numba as nb # type: ignore[import-not-found] +except ImportError: # pragma: no cover - performance-only optional dependency + + class _NumbaCompat: + """Fallback decorator shim when numba is not installed.""" + + @staticmethod + def njit(*args, **kwargs): # noqa: ARG004 + if args and callable(args[0]): + return args[0] + + def decorator(fn): + return fn + + return decorator + + nb = _NumbaCompat() + + +_COMPLEX_ZERO = 0j +_PI = np.pi +_E = np.e +_T_PHASE = np.exp(1j * _PI / 4) +_TDG_PHASE = np.exp(-1j * _PI / 4) + + +@nb.njit(cache=True, nogil=True, fastmath=False) +def _rz_phases(theta: float) -> tuple[complex, complex]: + """Parameterized Rz gate as diagonal phases.""" + half_theta = theta / 2 + return np.exp(-1j * half_theta), np.exp(1j * half_theta) + + +@nb.njit(cache=True, nogil=True, fastmath=False) +def ry(theta: float) -> np.ndarray: + """Parameterized Ry gate.""" + c, s = np.cos(theta / 2), np.sin(theta / 2) + mat = np.empty((2, 2), dtype=np.complex128) + mat[0, 0] = c + mat[0, 1] = -s + mat[1, 0] = s + mat[1, 1] = c + return mat + + +@nb.njit(cache=True, nogil=True, fastmath=False) +def rx(theta: float) -> np.ndarray: + """Parameterized Rx gate.""" + c, s = np.cos(theta / 2), np.sin(theta / 2) + mat = np.empty((2, 2), dtype=np.complex128) + mat[0, 0] = c + mat[0, 1] = -1j * s + mat[1, 0] = -1j * s + mat[1, 1] = c + return mat + + +def rz(theta: float) -> np.ndarray: + """Parameterized Rz gate.""" + phase0, phase1 = _rz_phases(theta) + return _diag_to_matrix(phase0, phase1) + + +@nb.njit(cache=True, nogil=True, fastmath=False) +def u3(theta: float, phi: float, lam: float) -> np.ndarray: + """Parameterized U3 gate (generic single-qubit rotation).""" + c, s = np.cos(theta / 2), np.sin(theta / 2) + exp_phi = np.exp(1j * phi) + exp_lam = np.exp(1j * lam) + exp_philam = np.exp(1j * (phi + lam)) + mat = np.empty((2, 2), dtype=np.complex128) + mat[0, 0] = c + mat[0, 1] = -exp_lam * s + mat[1, 0] = exp_phi * s + mat[1, 1] = exp_philam * c + return mat + + +@nb.njit(cache=True, nogil=True, fastmath=False) +def _diag_to_matrix(phase0: complex, phase1: complex) -> np.ndarray: + """Convert diagonal phases to a 2x2 matrix.""" + mat = np.zeros((2, 2), dtype=np.complex128) + mat[0, 0] = phase0 + mat[1, 1] = phase1 + return mat + + +@nb.njit(cache=True, nogil=True, fastmath=False) +def _is_diagonal_matrix(mat: np.ndarray, tol: float = 1e-10) -> bool: + """Return whether a 2x2 matrix is diagonal within tolerance.""" + return (np.abs(mat[0, 1]) < tol) and (np.abs(mat[1, 0]) < tol) + + +@nb.njit(cache=True, nogil=True, fastmath=False) +def _diagonal_phases_from_matrix(mat: np.ndarray) -> tuple[complex, complex]: + """Extract diagonal phases from a diagonal 2x2 matrix.""" + return mat[0, 0], mat[1, 1] + + +PARAMETERIZED_GATES = { + "rz": rz, + "ry": ry, + "rx": rx, + "u3": u3, +} + +_CONST_X = np.array([[0, 1], [1, 0]], dtype=np.complex128) +_CONST_Y = np.array([[0, -1j], [1j, 0]], dtype=np.complex128) +_CONST_Z = np.array([[1, 0], [0, -1]], dtype=np.complex128) +_CONST_H = np.array([[1, 1], [1, -1]], dtype=np.complex128) / np.sqrt(2) +_CONST_ID = np.eye(2, dtype=np.complex128) +_CONST_S = np.array([[1, 0], [0, 1j]], dtype=np.complex128) +_CONST_T = np.array([[1, 0], [0, _T_PHASE]], dtype=np.complex128) +_CONST_SDG = np.array([[1, 0], [0, -1j]], dtype=np.complex128) +_CONST_TDG = np.array([[1, 0], [0, _TDG_PHASE]], dtype=np.complex128) +_CONST_SX = 0.5 * np.array([[1 + 1j, 1 - 1j], [1 - 1j, 1 + 1j]], dtype=np.complex128) +_CONST_SXDG = 0.5 * np.array([[1 - 1j, 1 + 1j], [1 + 1j, 1 - 1j]], dtype=np.complex128) +_CONST_SWAP = np.array( + [[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]], dtype=np.complex128 +) + +NON_PARAMETERIZED_GATES: dict[str, np.ndarray] = { + "x": _CONST_X, + "y": _CONST_Y, + "z": _CONST_Z, + "h": _CONST_H, + "id": _CONST_ID, + "s": _CONST_S, + "t": _CONST_T, + "sdg": _CONST_SDG, + "tdg": _CONST_TDG, + "sx": _CONST_SX, + "sxdg": _CONST_SXDG, + "swap": _CONST_SWAP, +} + +# Fast-path single-qubit gate properties. Diagonal gates are kept as phase +# pairs while fusing so they can flush directly to the diagonal Cython kernel +# instead of materializing a 2x2 matrix. +NON_PARAMETERIZED_GATE_PROPS: dict[str, tuple[bool, complex, complex, np.ndarray | None]] = { + "x": (False, _COMPLEX_ZERO, _COMPLEX_ZERO, _CONST_X), + "y": (False, _COMPLEX_ZERO, _COMPLEX_ZERO, _CONST_Y), + "z": (True, 1.0 + _COMPLEX_ZERO, -1.0 + _COMPLEX_ZERO, None), + "h": (False, _COMPLEX_ZERO, _COMPLEX_ZERO, _CONST_H), + "id": (True, 1.0 + _COMPLEX_ZERO, 1.0 + _COMPLEX_ZERO, None), + "s": (True, 1.0 + _COMPLEX_ZERO, 1j, None), + "t": (True, 1.0 + _COMPLEX_ZERO, _T_PHASE, None), + "sdg": (True, 1.0 + _COMPLEX_ZERO, -1j, None), + "tdg": (True, 1.0 + _COMPLEX_ZERO, _TDG_PHASE, None), + "sx": (False, _COMPLEX_ZERO, _COMPLEX_ZERO, _CONST_SX), + "sxdg": (False, _COMPLEX_ZERO, _COMPLEX_ZERO, _CONST_SXDG), +} + +CONTROLLED_GATE_SUB_UNITARIES: dict[str, np.ndarray] = { + "cx": _CONST_X, + "cy": _CONST_Y, + "cz": _CONST_Z, +} + +# Pre-flatten non-parameterized gates for Cython (contiguous 1D arrays) +GATE_CACHE: dict[str, np.ndarray] = {} +for _name, _mat in NON_PARAMETERIZED_GATES.items(): + GATE_CACHE[_name] = np.ascontiguousarray(_mat.ravel()) +for _name, _mat in CONTROLLED_GATE_SUB_UNITARIES.items(): + GATE_CACHE[_name] = np.ascontiguousarray(_mat.ravel()) + +# Diagonal gate phases: gate_name -> (phase0, phase1) for diag(phase0, phase1) +DIAGONAL_PHASES: dict[str, tuple[complex, complex]] = { + "z": (1.0, -1.0), + "s": (1.0, 1j), + "t": (1.0, _T_PHASE), + "sdg": (1.0, -1j), + "tdg": (1.0, _TDG_PHASE), + "id": (1.0, 1.0), +} + +# Controlled diagonal phases: gate_name -> phase (applied when control=1 AND target=1) +CONTROLLED_DIAGONAL_PHASES: dict[str, complex] = { + "cz": -1.0, +} + +# Integer opcodes for preprocessed instructions +_OP_SINGLE = 0 +_OP_CONTROLLED = 1 +_OP_DIAGONAL = 2 +_OP_CTRL_DIAGONAL = 3 +_OP_TWO_QUBIT = 4 + + +def _try_eval_expression(expr) -> float | None: + """Try to evaluate an AST expression to a float value.""" + if isinstance(expr, (IntegerLiteral, FloatLiteral)): + return float(expr.value) + if isinstance(expr, Identifier): + if expr.name == "pi": + return _PI + if expr.name == "tau": + return 2 * _PI + if expr.name == "euler": + return _E + return None + if isinstance(expr, BinaryExpression): + lhs = _try_eval_expression(expr.lhs) + rhs = _try_eval_expression(expr.rhs) + if lhs is None or rhs is None: + return None + if expr.op is BinaryOperator["+"]: + return lhs + rhs + if expr.op is BinaryOperator["-"]: + return lhs - rhs + if expr.op is BinaryOperator["*"]: + return lhs * rhs + if expr.op is BinaryOperator["/"]: + return lhs / rhs + if expr.op is BinaryOperator["**"]: + return lhs**rhs + return None + if isinstance(expr, UnaryExpression): + operand = _try_eval_expression(expr.expression) + if operand is None: + return None + if expr.op is UnaryOperator["-"]: + return -operand + return None + return None + + +def _extract_params(statement): + """Extract float parameters from a QuantumGate's arguments.""" + params = [] + for arg in statement.arguments: + val = _try_eval_expression(arg) + if val is not None: + params.append(val) + return params + + +_CCX_ALIASES = frozenset({"ccx", "toffoli", "ccnot"}) + +# Two-qubit gates that are diagonal in the computational basis: pending +# diagonal one-qubit gates commute with these and need not be flushed. +_DIAGONAL_TWO_QUBIT_GATES = frozenset({"cz", "crz", "cp", "cphaseshift", "cu1"}) + +_PHASE_GATE_ALIASES = frozenset({"p", "phaseshift", "u1"}) +_CPHASE_GATE_ALIASES = frozenset({"cp", "cphaseshift", "cu1"}) + + +def _preprocess(program, num_qubits): + """Walk the unrolled AST once, producing packed numpy arrays with inline gate fusion. + + This keeps the existing no-cache semantics: every call preprocesses its + input program. The speedup comes from avoiding Python list growth, + preserving diagonal gates as phase pairs while fusing, and using optional + numba-compiled helpers for small matrix/phase constructors. + + Qubit operands are resolved through a ``(register name, local index) -> + global index`` map built from the ``QubitDeclaration`` statements of the + same AST, so multi-register programs address the correct amplitudes. + Statements the simulator cannot honor raise instead of being silently + dropped. + """ + statements = program._unrolled_ast.statements + if len(statements) == 0: + raise ValueError( + "Program has no unrolled statements. Call module.unroll() before simulating." + ) + + # Conservative bound: pending one-qubit flushes plus SWAP expansion can + # emit more instructions than source statements. The bound is deliberately + # simple and over-allocates a little to avoid dynamic list growth. + # ccx expands to 15 instructions, so those statements get an extra budget. + n_ccx = sum(1 for s in statements if isinstance(s, QuantumGate) and s.name.name in _CCX_ALIASES) + max_instructions = max(1, len(statements) * 4 + num_qubits + n_ccx * 15) + opcodes = np.empty(max_instructions, dtype=np.int32) + targets = np.empty(max_instructions, dtype=np.int32) + controls = np.empty(max_instructions, dtype=np.int32) + gate_params = np.empty(max_instructions * 4, dtype=np.complex128) + diag_phases = np.empty(max_instructions * 2, dtype=np.complex128) + two_qubit_offsets = np.empty(max_instructions, dtype=np.int32) + # The generic 4x4 two-qubit path is currently unused (swap decomposes to + # controlled-X; everything else maps to the specialized kernels), so its + # scratch buffer stays minimal. _write_two_qubit's slice assignment raises + # ValueError rather than silently overflowing if that ever changes. + two_qubit_gates = np.empty(16, dtype=np.complex128) + + # register name -> (global offset, size), built from the declarations in + # the same statement list we are about to walk. + reg_offsets: dict[str, tuple[int, int]] = {} + next_offset = 0 + + def _resolve_qubit(operand) -> int: + """Resolve a gate operand to its global qubit index.""" + if not isinstance(operand, IndexedIdentifier): + opname = getattr(operand, "name", None) + if isinstance(opname, str) and opname.startswith("$"): + raise NotImplementedError( + "Physical qubit operands (e.g. '$0') are not supported by the " + "statevector simulator." + ) + raise ValueError( + f"Unsupported qubit operand {operand!r}. Call module.unroll() before simulating." + ) + reg_name = operand.name.name + if reg_name not in reg_offsets: + raise ValueError(f"Gate references undeclared qubit register '{reg_name}'.") + offset, size = reg_offsets[reg_name] + first_index = operand.indices[0] + if not isinstance(first_index, list) or len(first_index) != 1: + raise ValueError( + f"Unsupported qubit index on register '{reg_name}'. " + "Call module.unroll() before simulating." + ) + index_node = first_index[0] + if not isinstance(index_node, IntegerLiteral): + raise ValueError( + f"Non-literal qubit index on register '{reg_name}'. " + "Call module.unroll() before simulating." + ) + local = index_node.value + if not 0 <= local < size: + raise ValueError( + f"Qubit index {local} out of range for register '{reg_name}' of size {size}." + ) + global_index = offset + local + if global_index >= num_qubits: + raise ValueError( + f"Resolved qubit index {global_index} exceeds the program's " + f"{num_qubits} qubits; the module's register bookkeeping is inconsistent." + ) + return global_index + + zero4 = np.zeros(4, dtype=np.complex128) + zero2 = np.zeros(2, dtype=np.complex128) + instruction_idx = 0 + tq_offset = 0 + + # target -> (is_diagonal, phase0, phase1, matrix) + pending: dict[int, tuple[bool, complex, complex, np.ndarray | None]] = {} + + def _write_single(target: int, mat: np.ndarray): + nonlocal instruction_idx + opcodes[instruction_idx] = _OP_SINGLE + targets[instruction_idx] = target + controls[instruction_idx] = -1 + gate_params[instruction_idx * 4 : instruction_idx * 4 + 4] = mat.ravel() + diag_phases[instruction_idx * 2 : instruction_idx * 2 + 2] = zero2 + two_qubit_offsets[instruction_idx] = -1 + instruction_idx += 1 + + def _write_diagonal(target: int, phase0: complex, phase1: complex): + nonlocal instruction_idx + opcodes[instruction_idx] = _OP_DIAGONAL + targets[instruction_idx] = target + controls[instruction_idx] = -1 + gate_params[instruction_idx * 4 : instruction_idx * 4 + 4] = zero4 + diag_phases[instruction_idx * 2 : instruction_idx * 2 + 2] = (phase0, phase1) + two_qubit_offsets[instruction_idx] = -1 + instruction_idx += 1 + + def _write_controlled(control: int, target: int, flat: np.ndarray): + nonlocal instruction_idx + opcodes[instruction_idx] = _OP_CONTROLLED + targets[instruction_idx] = target + controls[instruction_idx] = control + gate_params[instruction_idx * 4 : instruction_idx * 4 + 4] = flat + diag_phases[instruction_idx * 2 : instruction_idx * 2 + 2] = zero2 + two_qubit_offsets[instruction_idx] = -1 + instruction_idx += 1 + + def _write_controlled_diagonal(control: int, target: int, phase: complex): + nonlocal instruction_idx + opcodes[instruction_idx] = _OP_CTRL_DIAGONAL + targets[instruction_idx] = target + controls[instruction_idx] = control + gate_params[instruction_idx * 4 : instruction_idx * 4 + 4] = zero4 + diag_phases[instruction_idx * 2 : instruction_idx * 2 + 2] = (phase, 0j) + two_qubit_offsets[instruction_idx] = -1 + instruction_idx += 1 + + def _write_two_qubit(control: int, target: int, flat: np.ndarray): + nonlocal instruction_idx, tq_offset + opcodes[instruction_idx] = _OP_TWO_QUBIT + targets[instruction_idx] = target + controls[instruction_idx] = control + gate_params[instruction_idx * 4 : instruction_idx * 4 + 4] = zero4 + diag_phases[instruction_idx * 2 : instruction_idx * 2 + 2] = zero2 + two_qubit_offsets[instruction_idx] = tq_offset + two_qubit_gates[tq_offset : tq_offset + 16] = flat + tq_offset += 16 + instruction_idx += 1 + + def _write_ccx(ctrl_a: int, ctrl_b: int, tgt: int): + """Emit the standard exact 15-gate ccx decomposition (Nielsen & Chuang). + + H(t); CX(b,t); Tdg(t); CX(a,t); T(t); CX(b,t); Tdg(t); CX(a,t); + T(b); T(t); H(t); CX(a,b); T(a); Tdg(b); CX(a,b) + """ + cx_flat = GATE_CACHE["cx"] + t_p0, t_p1 = DIAGONAL_PHASES["t"] + tdg_p0, tdg_p1 = DIAGONAL_PHASES["tdg"] + _write_single(tgt, _CONST_H) + _write_controlled(ctrl_b, tgt, cx_flat) + _write_diagonal(tgt, tdg_p0, tdg_p1) + _write_controlled(ctrl_a, tgt, cx_flat) + _write_diagonal(tgt, t_p0, t_p1) + _write_controlled(ctrl_b, tgt, cx_flat) + _write_diagonal(tgt, tdg_p0, tdg_p1) + _write_controlled(ctrl_a, tgt, cx_flat) + _write_diagonal(ctrl_b, t_p0, t_p1) + _write_diagonal(tgt, t_p0, t_p1) + _write_single(tgt, _CONST_H) + _write_controlled(ctrl_a, ctrl_b, cx_flat) + _write_diagonal(ctrl_a, t_p0, t_p1) + _write_diagonal(ctrl_b, tdg_p0, tdg_p1) + _write_controlled(ctrl_a, ctrl_b, cx_flat) + + def _flush_pending(target: int): + if target not in pending: + return + is_diagonal, phase0, phase1, mat = pending.pop(target) + + if is_diagonal and (np.abs(phase0 - 1.0) < 1e-10) and (np.abs(phase1 - 1.0) < 1e-10): + return + if not is_diagonal and mat is not None and _is_diagonal_matrix(mat): + is_diagonal = True + phase0, phase1 = _diagonal_phases_from_matrix(mat) + mat = None + + if is_diagonal: + _write_diagonal(target, phase0, phase1) + else: + if mat is None: + mat = _diag_to_matrix(phase0, phase1) + _write_single(target, mat) + + def _flush_all(): + for target in list(pending.keys()): + _flush_pending(target) + + def _accumulate( + target: int, + new_is_diagonal: bool, + new_phase0: complex, + new_phase1: complex, + new_mat: np.ndarray | None, + ): + if target not in pending: + pending[target] = (new_is_diagonal, new_phase0, new_phase1, new_mat) + return + + current_is_diagonal, current_phase0, current_phase1, current_mat = pending[target] + if current_is_diagonal and new_is_diagonal: + pending[target] = ( + True, + new_phase0 * current_phase0, + new_phase1 * current_phase1, + None, + ) + return + + current_matrix = ( + _diag_to_matrix(current_phase0, current_phase1) if current_is_diagonal else current_mat + ) + new_matrix = _diag_to_matrix(new_phase0, new_phase1) if new_is_diagonal else new_mat + if current_matrix is None or new_matrix is None: + raise ValueError("Invalid pending gate fusion state.") + fused = new_matrix @ current_matrix + if _is_diagonal_matrix(fused): + phase0, phase1 = _diagonal_phases_from_matrix(fused) + pending[target] = (True, phase0, phase1, None) + else: + pending[target] = (False, _COMPLEX_ZERO, _COMPLEX_ZERO, fused) + + measured = False + for statement in statements: + if isinstance(statement, QubitDeclaration): + decl_name = statement.qubit.name + decl_size = statement.size.value if statement.size is not None else 1 + reg_offsets[decl_name] = (next_offset, decl_size) + next_offset += decl_size + continue + if isinstance(statement, (Include, ClassicalDeclaration, QuantumBarrier)): + # No effect on the statevector. + continue + if isinstance(statement, QuantumMeasurementStatement): + # Terminal measurement is honored by the sampling step in run(). + # Anything acting on the state afterwards makes it mid-circuit, + # which the simulator cannot honor (checked below). + measured = True + continue + if isinstance(statement, QuantumReset): + raise NotImplementedError("'reset' is not supported by the statevector simulator.") + if isinstance(statement, BranchingStatement): + raise NotImplementedError( + "Classical control flow ('if') is not supported by the statevector simulator." + ) + if not isinstance(statement, QuantumGate): + raise NotImplementedError( + f"Statement type '{type(statement).__name__}' is not supported " + "by the statevector simulator." + ) + if measured: + raise NotImplementedError( + "Mid-circuit measurement is not supported by the statevector " + "simulator: a gate acts on the state after a measurement." + ) + + gate_name = statement.name.name + params = _extract_params(statement) + + if len(statement.qubits) == 1: + target = _resolve_qubit(statement.qubits[0]) + + if gate_name in NON_PARAMETERIZED_GATE_PROPS: + is_diagonal, phase0, phase1, mat = NON_PARAMETERIZED_GATE_PROPS[gate_name] + _accumulate(target, is_diagonal, phase0, phase1, mat) + elif gate_name == "rz" and len(params) == 1: + phase0, phase1 = _rz_phases(params[0]) + _accumulate(target, True, phase0, phase1, None) + elif gate_name in _PHASE_GATE_ALIASES and len(params) == 1: + # p(theta) == diag(1, e^{i*theta}) + _accumulate(target, True, 1.0 + _COMPLEX_ZERO, np.exp(1j * params[0]), None) + elif gate_name in PARAMETERIZED_GATES: + gate_fn = PARAMETERIZED_GATES[gate_name] + required_params = 3 if gate_name == "u3" else 1 + if len(params) != required_params: + raise ValueError(f"Gate {gate_name} requires {required_params} parameter(s).") + mat = gate_fn(*params) + _accumulate(target, False, _COMPLEX_ZERO, _COMPLEX_ZERO, mat) + else: + raise ValueError( + f"Gate '{gate_name}' is not supported by the statevector simulator. " + "If the program has not been unrolled, call module.unroll() first." + ) + + elif len(statement.qubits) == 2: + control = _resolve_qubit(statement.qubits[0]) + target = _resolve_qubit(statement.qubits[1]) + if control == target: + raise ValueError( + f"Gate '{gate_name}' has identical control and target " + f"qubit (global index {target})." + ) + + # Flush pending gates only when they do not commute with the + # incoming two-qubit gate. Diagonal one-qubit gates commute with + # diagonal two-qubit gates (CZ/CRZ/CP), so they can remain fused + # until a later non-diagonal interaction or the final flush. + for qubit in (control, target): + if qubit not in pending: + continue + is_diagonal, _, _, _ = pending[qubit] + if (not is_diagonal) or gate_name not in _DIAGONAL_TWO_QUBIT_GATES: + _flush_pending(qubit) + + if gate_name == "swap": + # Decompose SWAP into three controlled-X instructions. This + # stays on the optimized controlled-gate kernel and avoids the + # generic 4x4 two-qubit path for this common gate. + flat = GATE_CACHE["cx"] + _write_controlled(control, target, flat) + _write_controlled(target, control, flat) + _write_controlled(control, target, flat) + elif gate_name in CONTROLLED_DIAGONAL_PHASES: + phase = CONTROLLED_DIAGONAL_PHASES[gate_name] + _write_controlled_diagonal(control, target, phase) + elif gate_name in _CPHASE_GATE_ALIASES and len(params) == 1: + # cp(theta) phases only the |11> amplitude: exactly the + # controlled-diagonal kernel. + _write_controlled_diagonal(control, target, np.exp(1j * params[0])) + elif gate_name == "crz" and len(params) == 1: + # CRZ applies Rz(theta) = diag(e^{-i*theta/2}, e^{i*theta/2}) to + # the target when the control is set, so BOTH target=0 and + # target=1 pick up a phase. A controlled-diagonal instruction + # can only phase the |control=1, target=1> amplitude, which would + # implement a controlled-phase gate instead. Route it through the + # general controlled-gate kernel with the full (diagonal) Rz. + theta = params[0] + phase0, phase1 = _rz_phases(theta) + flat = np.array([phase0, _COMPLEX_ZERO, _COMPLEX_ZERO, phase1], dtype=np.complex128) + _write_controlled(control, target, flat) + elif gate_name in CONTROLLED_GATE_SUB_UNITARIES: + flat = GATE_CACHE[gate_name] + _write_controlled(control, target, flat) + else: + raise ValueError( + f"Gate '{gate_name}' is not supported by the statevector simulator. " + "If the program has not been unrolled, call module.unroll() first." + ) + + elif len(statement.qubits) == 3 and gate_name in _CCX_ALIASES: + ctrl_a = _resolve_qubit(statement.qubits[0]) + ctrl_b = _resolve_qubit(statement.qubits[1]) + tgt = _resolve_qubit(statement.qubits[2]) + if len({ctrl_a, ctrl_b, tgt}) != 3: + raise ValueError( + f"Gate '{gate_name}' requires three distinct qubits, got " + f"global indices ({ctrl_a}, {ctrl_b}, {tgt})." + ) + # The decomposition mixes diagonal and non-diagonal gates on all + # three qubits, so flush anything pending on them first. + for qubit in (ctrl_a, ctrl_b, tgt): + _flush_pending(qubit) + _write_ccx(ctrl_a, ctrl_b, tgt) + + else: + raise ValueError( + f"Gate '{gate_name}' acting on {len(statement.qubits)} qubits " + "is not supported by the statevector simulator." + ) + + # Flush remaining pending single-qubit gates + _flush_all() + + n = instruction_idx + if n == 0: + return n, None, None, None, None, None, None, None + + if tq_offset > 0: + two_qubit_gates_trimmed = np.ascontiguousarray(two_qubit_gates[:tq_offset]) + else: + two_qubit_gates_trimmed = np.zeros(1, dtype=np.complex128) + + return ( + n, + opcodes[:n], + targets[:n], + controls[:n], + np.ascontiguousarray(gate_params[: n * 4]), + np.ascontiguousarray(diag_phases[: n * 2]), + two_qubit_offsets[:n], + two_qubit_gates_trimmed, + ) + + +@dataclass(frozen=True) +class SimulatorResult: + """Class to store the result of a statevector simulation. + + Conventions (matching qiskit): ``probabilities`` and ``final_statevector`` + are indexed little-endian, i.e. bit ``k`` of the array index is the state + of qubit ``k``. ``measurement_counts`` keys are the binary rendering of + that same index, so qubit 0 is the RIGHTMOST character and + ``probabilities[int(key, 2)]`` is the probability of outcome ``key``. + """ + + probabilities: np.ndarray + measurement_counts: Counter[str] + final_statevector: np.ndarray + + +# Statevector memory ceiling: 2**n complex128 amplitudes = 16 * 2**n bytes. +# The default of 30 qubits caps the allocation at 16 GiB; without a guard the +# lazily-committed allocation dies with SIGBUS/OOM mid-kernel instead of +# raising. Override with PYQASM_SIM_MAX_QUBITS for larger machines. +_DEFAULT_MAX_QUBITS = 30 + + +def _max_sim_qubits() -> int: + """Resolve the simulator qubit ceiling from the environment.""" + raw = os.environ.get("PYQASM_SIM_MAX_QUBITS", "") + try: + value = int(raw) + except ValueError: + return _DEFAULT_MAX_QUBITS + return value if value > 0 else _DEFAULT_MAX_QUBITS + + +class Simulator: + """ + Statevector simulator + + """ + + def __init__(self, seed: None | int = None): + """ + Initialize the statevector simulator. + + Parameters: + seed (None | int): A seed to initialize the `np.random.BitGenerator`. + If None, then fresh, unpredictable entropy will be pulled from + the OS using `np.random.SeedSequence`. + """ + self._rng = np.random.default_rng(seed) + + def run(self, program: QasmModule | str, shots: int = 1) -> SimulatorResult: + """Run the statevector simulator. + + For QasmModule input, the module must already be unrolled (call + module.unroll(), and optionally module.remove_idle_qubits(), + beforehand); an un-unrolled module raises rather than silently + simulating the raw AST. + + For string input, loads/unrolls/removes idle qubits automatically. + + Args: + program (QasmModule | str): The program to simulate. + shots (int): The number of shots to simulate. Defaults to 1. + + Returns: + SimulatorResult: The result of the simulation. See + SimulatorResult for the bit-ordering conventions. + + Raises: + ValueError: If shots is less than 0, if the program contains an + unsupported gate or invalid syntax, if a QasmModule input has + not been unrolled, or if the program exceeds the qubit ceiling + (default 30; override with PYQASM_SIM_MAX_QUBITS). + NotImplementedError: If the program uses reset, mid-circuit + measurement, classical control flow, or physical qubits. + """ + if shots < 0: + raise ValueError("Shots must be greater than or equal to 0.") + + if isinstance(program, str): + program = loads(program) + program.unroll() + program.remove_idle_qubits() + elif len(program._unrolled_ast.statements) == 0: + raise ValueError( + "QasmModule has not been unrolled. Call module.unroll() before Simulator.run()." + ) + + num_qubits = program.num_qubits + if num_qubits == 0: + trivial = np.ones(1, dtype=np.complex128) + counts: Counter[str] = Counter({"": shots}) if shots > 0 else Counter() + return SimulatorResult(np.abs(trivial) ** 2, counts, trivial) + + max_qubits = _max_sim_qubits() + if num_qubits > max_qubits: + required_gib = (16 * 2**num_qubits) / 2**30 + raise ValueError( + f"Program has {num_qubits} qubits; the statevector alone would " + f"require {required_gib:.0f} GiB. The simulator's ceiling is " + f"{max_qubits} qubits (override with PYQASM_SIM_MAX_QUBITS)." + ) + + n, opcodes, targets, controls, gate_params, diag_phases, tq_offsets, tq_gates = _preprocess( + program, num_qubits + ) + + sv = np.zeros(2**num_qubits, dtype=np.complex128) + sv[0] = 1.0 + + if n > 0: + apply_circuit( + sv, + num_qubits, + opcodes, + targets, + controls, + gate_params, + diag_phases, + tq_offsets, + tq_gates, + n, + ) + + probabilities = np.abs(sv) ** 2 + samples = self._rng.choice(len(probabilities), size=shots, p=probabilities) + # Little-endian, qiskit-compatible keys: bit k of the index is qubit k, + # so qubit 0 is the rightmost character (see SimulatorResult). + counts = Counter(format(s, f"0{num_qubits}b") for s in samples) + + return SimulatorResult(probabilities, counts, sv) diff --git a/tests/test_perf_regression.py b/tests/test_perf_regression.py new file mode 100644 index 00000000..ab56caf7 --- /dev/null +++ b/tests/test_perf_regression.py @@ -0,0 +1,186 @@ +# Copyright 2025 qBraid +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# pylint: disable=redefined-outer-name + +""" +Performance regression tests for the PyQASM statevector simulator. + +These tests enforce time budgets to catch performance regressions. Each test +runs the simulator multiple times and checks that the median execution time +stays within a generous upper bound (3x the baseline measured on Apple Silicon). + +Run with: pytest tests/test_perf_regression.py -v +Skip in CI: pytest -m "not benchmark" + +Time budgets are intentionally loose (3x baseline) to avoid flakiness on +different hardware. The goal is catching order-of-magnitude regressions, +not enforcing tight timings. + +Baselines (Apple Silicon M-series, March 2026): + Random 16q/200d: ~24 ms + Random 20q/300d: ~330 ms + Random 22q/400d: ~1625 ms + QFT 10q: ~27 ms + QFT 14q: ~58 ms +""" + +import random +import time + +import numpy as np +import pytest + +from pyqasm import loads as pyqasm_loads +from pyqasm.simulator.statevector import Simulator + +SINGLE_QUBIT_GATES = ["h", "x", "y", "z", "s", "t", "rx(1.0)", "ry(0.5)", "rz(0.3)"] +TWO_QUBIT_GATES = ["cx", "cy", "cz", "swap"] + +N_REPEATS = 5 +BUDGET_MULTIPLIER = 3 # allow 3x baseline for portability across hardware + + +def generate_random_qasm(num_qubits: int, depth: int, seed: int = 42) -> str: + rng = random.Random(seed) + lines = [ + "OPENQASM 3;", + 'include "stdgates.inc";', + f"qubit[{num_qubits}] q;", + ] + for _ in range(depth): + if num_qubits >= 2 and rng.random() < 0.4: + gate = rng.choice(TWO_QUBIT_GATES) + q0, q1 = rng.sample(range(num_qubits), 2) + lines.append(f"{gate} q[{q0}], q[{q1}];") + else: + gate = rng.choice(SINGLE_QUBIT_GATES) + q = rng.randint(0, num_qubits - 1) + lines.append(f"{gate} q[{q}];") + return "\n".join(lines) + + +def generate_qft_qasm(num_qubits: int) -> str: + lines = [ + "OPENQASM 3;", + 'include "stdgates.inc";', + f"qubit[{num_qubits}] q;", + ] + for i in range(num_qubits): + lines.append(f"h q[{i}];") + for j in range(i + 1, num_qubits): + k = j - i + angle = f"pi/{2**k}" + lines.append(f"crz({angle}) q[{j}], q[{i}];") + for i in range(num_qubits // 2): + lines.append(f"swap q[{i}], q[{num_qubits - 1 - i}];") + return "\n".join(lines) + + +def median_sim_time(qasm: str, sim: Simulator, n_repeats: int = N_REPEATS) -> float: + """Return the median wall-clock time (seconds) of sim.run() over n_repeats.""" + module = pyqasm_loads(qasm) + module.unroll() + module.remove_idle_qubits() + times = [] + for _ in range(n_repeats): + start = time.perf_counter() + sim.run(module, shots=0) + times.append(time.perf_counter() - start) + return float(np.median(times)) + + +@pytest.fixture +def sim(): + return Simulator(seed=0) + + +# --------------------------------------------------------------------------- +# Random circuit performance tests +# --------------------------------------------------------------------------- + + +@pytest.mark.benchmark +@pytest.mark.parametrize( + "num_qubits, depth, baseline_ms", + [ + (16, 200, 24), + (20, 300, 330), + ], + ids=["random-16q", "random-20q"], +) +def test_random_circuit_perf(sim, num_qubits, depth, baseline_ms): + """Median simulation time must stay within budget for random circuits.""" + qasm = generate_random_qasm(num_qubits, depth) + elapsed = median_sim_time(qasm, sim) + budget_ms = baseline_ms * BUDGET_MULTIPLIER + elapsed_ms = elapsed * 1000 + assert elapsed_ms < budget_ms, ( + f"Random {num_qubits}q/{depth}d: {elapsed_ms:.1f} ms exceeded " + f"budget {budget_ms:.0f} ms (baseline {baseline_ms} ms x{BUDGET_MULTIPLIER})" + ) + + +# --------------------------------------------------------------------------- +# QFT circuit performance tests +# --------------------------------------------------------------------------- + + +@pytest.mark.benchmark +@pytest.mark.parametrize( + "num_qubits, baseline_ms", + [ + (10, 27), + (14, 58), + ], + ids=["qft-10q", "qft-14q"], +) +def test_qft_circuit_perf(sim, num_qubits, baseline_ms): + """Median simulation time must stay within budget for QFT circuits.""" + qasm = generate_qft_qasm(num_qubits) + elapsed = median_sim_time(qasm, sim) + budget_ms = baseline_ms * BUDGET_MULTIPLIER + elapsed_ms = elapsed * 1000 + assert elapsed_ms < budget_ms, ( + f"QFT {num_qubits}q: {elapsed_ms:.1f} ms exceeded " + f"budget {budget_ms:.0f} ms (baseline {baseline_ms} ms x{BUDGET_MULTIPLIER})" + ) + + +# --------------------------------------------------------------------------- +# Scaling sanity check: ensure O(2^n) per gate, not worse +# --------------------------------------------------------------------------- + + +@pytest.mark.benchmark +def test_scaling_not_superexponential(sim): + """Verify that doubling qubits roughly doubles per-gate time (not worse). + + Compares 14q vs 16q random circuits with the same gate count. The ratio + of per-gate times should be near 4x (since state size quadruples with +2 + qubits). We allow up to 8x to account for cache effects and noise. + """ + depth = 150 + qasm_14 = generate_random_qasm(14, depth) + qasm_16 = generate_random_qasm(16, depth) + + t14 = median_sim_time(qasm_14, sim) + t16 = median_sim_time(qasm_16, sim) + + ratio = t16 / t14 + # +2 qubits => 4x state size => expect ~4x time. Allow up to 8x. + assert ratio < 8.0, ( + f"Scaling ratio 16q/14q = {ratio:.1f}x, expected <8x " + f"(14q={t14*1000:.1f}ms, 16q={t16*1000:.1f}ms)" + ) diff --git a/tests/test_sv_sim.py b/tests/test_sv_sim.py new file mode 100644 index 00000000..700d3478 --- /dev/null +++ b/tests/test_sv_sim.py @@ -0,0 +1,438 @@ +# Copyright 2025 qBraid +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +# pylint: disable=redefined-outer-name + +""" +Module containing tests for the PyQASM statevector simulator. + +""" + +import numpy as np +import pytest + +# Qiskit is only required for these comparison tests and is not installable in +# every CI environment (e.g. the manylinux wheel-build containers), so skip the +# whole module gracefully when it is unavailable. +pytest.importorskip("qiskit") +pytest.importorskip("qiskit_aer") + +# pylint: disable=wrong-import-position +from qiskit import transpile # noqa: E402 +from qiskit.qasm3 import loads # noqa: E402 +from qiskit_aer import AerSimulator # noqa: E402 + +from pyqasm import loads as pyqasm_loads # noqa: E402 +from pyqasm.simulator.statevector import Simulator # noqa: E402 + +# pylint: enable=wrong-import-position + + +@pytest.fixture +def aer_simulator(): + """Fixture to create and return a Qiskit AerSimulator.""" + return AerSimulator(method="statevector") + + +@pytest.fixture +def pyqasm_simulator(): + """Fixture to create and return a pyqasm simulator.""" + return Simulator(seed=22) + + +@pytest.mark.parametrize( + "qasm, description", + [ + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[1] q; + h q[0]; + """, + "Hadamard gate", + ), + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[1] q; + rx(pi) q[0]; + """, + "RX(pi) gate", + ), + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[1] q; + ry(pi) q[0]; + """, + "RY(pi) gate", + ), + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[1] q; + rz(pi) q[0]; + """, + "RZ(pi) gate", + ), + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[1] q; + s q[0]; + """, + "S gate", + ), + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[1] q; + t q[0]; + """, + "T gate", + ), + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[1] q; + sdg q[0]; + """, + "Sdg gate", + ), + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[1] q; + tdg q[0]; + """, + "Tdg gate", + ), + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[1] q; + id q[0]; + """, + "Identity gate", + ), + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[1] q; + u3(1.0, 0.5, 0.3) q[0]; + """, + "U3 gate (global phase)", + ), + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[1] q; + h q[0]; + rx(pi) q[0]; + ry(pi) q[0]; + rz(pi) q[0]; + s q[0]; + t q[0]; + sdg q[0]; + tdg q[0]; + id q[0]; + """, + "Combined gates: H, RX(pi), RY(pi), RZ(pi), S, T, Sdg, Tdg, ID", + ), + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[2] q; + cx q[0], q[1]; + """, + "CX gate", + ), + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[2] q; + x q[0]; + x q[1]; + cy q[0], q[1]; + """, + "CY on |11>", + ), + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[2] q; + x q[0]; + x q[1]; + cz q[0], q[1]; + """, + "CZ on |11>", + ), + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[2] q; + x q[1]; + swap q[0], q[1]; + """, + "SWAP on |01>", + ), + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[3] q; + x q[0]; + id q[1]; + id q[2]; + """, + "X on qubit 0 of 3-qubit system", + ), + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[2] q; + cy q[0], q[1]; + """, + "CY gate", + ), + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[2] q; + cz q[0], q[1]; + """, + "CZ gate", + ), + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[2] q; + swap q[0], q[1]; + """, + "Swap gate", + ), + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[2] q; + h q[0]; + crz(pi/4) q[0], q[1]; + """, + "CRZ gate", + ), + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[2] q; + h q[0]; + h q[1]; + """, + "Two H Gates", + ), + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[2] q; + h q[0]; + cx q[0], q[1]; + """, + "Bell state", + ), + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[3] q; + h q[0]; + cx q[0], q[1]; + cx q[1], q[2]; + """, + "GHZ state", + ), + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[2] q; + h q[0]; + swap q[0], q[1]; + cx q[0], q[1]; + """, + "Bell state with swap", + ), + ( + """ + OPENQASM 3; + include "stdgates.inc"; + qubit[2] q; + h q[0]; + swap q[1], q[0]; + cx q[0], q[1]; + """, + "Bell state with swap rev qubits", + ), + ], +) +def test_simulator_sv_results(qasm, description, aer_simulator, pyqasm_simulator): + """Test PyQASM simulator by comparing SV results against the Qiskit Aer simulator.""" + circuit = loads(qasm) + circuit.save_statevector() + compiled_circuit = transpile(circuit, aer_simulator, optimization_level=0) + result = aer_simulator.run(compiled_circuit).result() + sv_qiskit = result.get_statevector(compiled_circuit) + sv_expected = np.asarray(sv_qiskit) + + result = pyqasm_simulator.run(qasm, shots=1000) + + sv_actual = result.final_statevector + # Compare up to global phase unconditionally: equivalent states may differ + # by an overall e^{i*phi}, which is physically irrelevant. The unit-modulus + # assertion still catches any non-phase deviation in magnitude. + idx = int(np.argmax(np.abs(sv_expected))) + phase = sv_actual[idx] / sv_expected[idx] + assert np.isclose(abs(phase), 1.0), ( + f"Test failed for: {description}\n" + f"Statevectors differ by more than a global phase\n" + f"PyQASM simulator statevector: {result.final_statevector}\n" + f"Expected statevector: {sv_expected}" + ) + sv_actual = sv_actual / phase + + assert np.allclose(sv_actual, sv_expected), ( + f"Test failed for: {description}\n" + f"PyQASM simulator statevector: {result.final_statevector}\n" + f"Expected statevector: {sv_expected}" + ) + + +def test_simulator_large_circuit(aer_simulator, pyqasm_simulator): + """Test PyQASM simulator on an 8-qubit circuit.""" + qasm = """ + OPENQASM 3; + include "stdgates.inc"; + qubit[8] q; + h q[0]; + cx q[0], q[1]; + cx q[1], q[2]; + cx q[2], q[3]; + h q[4]; + cx q[4], q[5]; + cx q[5], q[6]; + cx q[6], q[7]; + swap q[3], q[4]; + """ + circuit = loads(qasm) + circuit.save_statevector() + compiled_circuit = transpile(circuit, aer_simulator, optimization_level=0) + result_qiskit = aer_simulator.run(compiled_circuit).result() + sv_expected = np.asarray(result_qiskit.get_statevector(compiled_circuit)) + + result = pyqasm_simulator.run(qasm, shots=0) + + assert np.allclose(result.final_statevector, sv_expected), ( + f"Test failed for 8-qubit circuit\n" + f"PyQASM: {result.final_statevector}\n" + f"Expected: {sv_expected}" + ) + + +def test_simulator_ccx_matches_qiskit(aer_simulator, pyqasm_simulator): + """ccx (previously silently dropped) must match qiskit exactly.""" + qasm = """ + OPENQASM 3; + include "stdgates.inc"; + qubit[3] q; + ry(0.7) q[0]; + ry(1.1) q[1]; + ry(0.4) q[2]; + ccx q[0], q[1], q[2]; + ccx q[2], q[0], q[1]; + """ + circuit = loads(qasm) + circuit.save_statevector() + compiled = transpile(circuit, aer_simulator, optimization_level=0) + sv_expected = np.asarray(aer_simulator.run(compiled).result().get_statevector(compiled)) + + result = pyqasm_simulator.run(qasm, shots=0) + assert np.allclose(result.final_statevector, sv_expected) + + +def test_simulator_multi_register_matches_qiskit(aer_simulator, pyqasm_simulator): + """Multi-register programs (previously mapped to the wrong qubits) must + match qiskit's declaration-order qubit layout.""" + qasm = """ + OPENQASM 3; + include "stdgates.inc"; + qubit[2] a; + qubit[2] b; + h a[0]; + cx a[0], b[1]; + x b[0]; + cz a[1], b[0]; + """ + circuit = loads(qasm) + circuit.save_statevector() + compiled = transpile(circuit, aer_simulator, optimization_level=0) + sv_expected = np.asarray(aer_simulator.run(compiled).result().get_statevector(compiled)) + + module = pyqasm_loads(qasm) + module.unroll() # no remove_idle_qubits: keep the full 4-qubit layout + result = pyqasm_simulator.run(module, shots=0) + assert np.allclose(result.final_statevector, sv_expected) + + +def test_simulator_counts_convention_matches_qiskit(aer_simulator, pyqasm_simulator): + """Counts keys use qiskit's convention: qubit 0 is the rightmost character.""" + qasm_measured = """ + OPENQASM 3; + include "stdgates.inc"; + qubit[3] q; + bit[3] c; + x q[2]; + c = measure q; + """ + circuit = loads(qasm_measured) + compiled = transpile(circuit, aer_simulator, optimization_level=0) + qiskit_counts = aer_simulator.run(compiled, shots=10).result().get_counts() + + module = pyqasm_loads( + 'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[3] q;\nid q[0];\nid q[1];\nx q[2];\n' + ) + module.unroll() + result = pyqasm_simulator.run(module, shots=10) + + assert list(qiskit_counts) == ["100"] + assert dict(result.measurement_counts) == {"100": 10} + assert result.probabilities[int("100", 2)] == pytest.approx(1.0) diff --git a/tests/test_sv_sim_core.py b/tests/test_sv_sim_core.py new file mode 100644 index 00000000..34579b68 --- /dev/null +++ b/tests/test_sv_sim_core.py @@ -0,0 +1,703 @@ +# Copyright 2025 qBraid +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +""" +Qiskit-free tests for the PyQASM statevector simulator. + +``test_sv_sim.py`` compares against Qiskit and is skipped wherever Qiskit is not +installed (e.g. the wheel-build containers). The tests here have no third-party +dependency so they run everywhere, and they focus on correctness of the +fast-path gate handling, fusion, and error behaviour rather than raw numerics. +""" + +import numpy as np +import pytest +from openqasm3.ast import ( + BinaryExpression, + BinaryOperator, + BooleanLiteral, + FloatLiteral, + Identifier, + IntegerLiteral, + UnaryExpression, + UnaryOperator, +) + +from pyqasm import loads +from pyqasm.simulator.statevector import ( + Simulator, + _try_eval_expression, + rz, +) + +# --- tiny reference statevector simulator (little-endian: qubit i is bit i) --- + +_I = np.eye(2, dtype=complex) +_REF_GATES = { + "x": np.array([[0, 1], [1, 0]], dtype=complex), + "y": np.array([[0, -1j], [1j, 0]], dtype=complex), + "z": np.array([[1, 0], [0, -1]], dtype=complex), + "h": np.array([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2), + "s": np.array([[1, 0], [0, 1j]], dtype=complex), + "sdg": np.array([[1, 0], [0, -1j]], dtype=complex), + "t": np.array([[1, 0], [0, np.exp(1j * np.pi / 4)]], dtype=complex), + "tdg": np.array([[1, 0], [0, np.exp(-1j * np.pi / 4)]], dtype=complex), + "id": _I, +} + + +def _embed(mat, target, num_qubits): + op = np.array([[1]], dtype=complex) + for qubit in range(num_qubits): + op = np.kron(mat if qubit == target else _I, op) + return op + + +def _embed_controlled(mat, control, target, num_qubits): + p0 = np.array([[1, 0], [0, 0]], dtype=complex) + p1 = np.array([[0, 0], [0, 1]], dtype=complex) + op0 = op1 = np.array([[1]], dtype=complex) + for qubit in range(num_qubits): + if qubit == control: + op0, op1 = np.kron(p0, op0), np.kron(p1, op1) + elif qubit == target: + op0, op1 = np.kron(_I, op0), np.kron(mat, op1) + else: + op0, op1 = np.kron(_I, op0), np.kron(_I, op1) + return op0 + op1 + + +def reference_statevector(ops, num_qubits): + """ops: list of (name, [qubits], [params]). Returns the final statevector.""" + sv = np.zeros(2**num_qubits, dtype=complex) + sv[0] = 1.0 + for name, qubits, params in ops: + if name in _REF_GATES: + sv = _embed(_REF_GATES[name], qubits[0], num_qubits) @ sv + elif name == "rz": + mat = np.array([[np.exp(-1j * params[0] / 2), 0], [0, np.exp(1j * params[0] / 2)]]) + sv = _embed(mat, qubits[0], num_qubits) @ sv + elif name in ("cx", "cy", "cz"): + sv = _embed_controlled(_REF_GATES[name[1:]], qubits[0], qubits[1], num_qubits) @ sv + elif name == "crz": + mat = np.array([[np.exp(-1j * params[0] / 2), 0], [0, np.exp(1j * params[0] / 2)]]) + sv = _embed_controlled(mat, qubits[0], qubits[1], num_qubits) @ sv + else: + raise ValueError(f"reference simulator has no gate {name!r}") + return sv + + +def assert_sv_close(actual, expected, atol=1e-9): + """Compare statevectors up to an irrelevant global phase.""" + idx = int(np.argmax(np.abs(expected))) + phase = actual[idx] / expected[idx] + assert abs(abs(phase) - 1) < atol, "statevectors differ by more than a global phase" + assert np.allclose(actual, phase * expected, atol=atol) + + +def run_sv(qasm, external_gates=None): + module = loads(qasm) + module.unroll(external_gates=external_gates) + return Simulator(seed=0).run(module, shots=0).final_statevector + + +# -------------------------------------------------------------------------- +# crz fast path regression: a controlled-Rz phases BOTH target=0 and target=1 +# when the control is set. The fast path must equal pyqasm's own (qiskit- +# verified) full decomposition; the earlier implementation only phased the +# |control=1, target=1> amplitude, which is a controlled-phase gate. +# -------------------------------------------------------------------------- + + +@pytest.mark.parametrize("theta", [0.0, 0.7, np.pi / 2, np.pi, -1.3, 2 * np.pi]) +@pytest.mark.parametrize("control, target", [(0, 1), (1, 0)]) +def test_crz_fast_path_matches_decomposition(theta, control, target): + qasm = ( + "OPENQASM 3;\n" + 'include "stdgates.inc";\n' + "qubit[2] q;\n" + "h q[0];\n" + "h q[1];\n" + f"crz({theta}) q[{control}], q[{target}];\n" + ) + fast = run_sv(qasm, external_gates=["crz"]) # hits the crz fast path + decomposed = run_sv(qasm) # full rz/rx/cx decomposition (oracle) + assert_sv_close(fast, decomposed) + + +def test_crz_phases_control_one_target_zero(): + """Directly pin the |control=1, target=0> phase the bug used to drop.""" + theta = 1.0 + qasm = ( + "OPENQASM 3;\n" + 'include "stdgates.inc";\n' + "qubit[2] q;\n" + "x q[0];\n" # control = 1, target = 0 -> basis index 1 + f"crz({theta}) q[0], q[1];\n" + ) + sv = run_sv(qasm, external_gates=["crz"]) + expected = reference_statevector([("x", [0], []), ("crz", [0, 1], [theta])], 2) + assert_sv_close(sv, expected) + # the amplitude at |control=1,target=0> (index 1) must carry e^{-i*theta/2} + assert np.isclose(sv[1], np.exp(-1j * theta / 2)) + + +# -------------------------------------------------------------------------- +# Known-action checks for the directly-supported fast-path gates. +# -------------------------------------------------------------------------- + + +@pytest.mark.parametrize("gate", ["x", "y", "z", "h", "s", "sdg", "t", "tdg", "id"]) +def test_single_qubit_gate_action(gate): + qasm = f'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[1] q;\nx q[0];\n{gate} q[0];\n' + sv = Simulator(seed=0).run(loads_unrolled(qasm), shots=0).final_statevector + expected = reference_statevector([("x", [0], []), (gate, [0], [])], 1) + assert_sv_close(sv, expected) + + +def test_cz_swap_cx_action(): + # cz on |11> -> -|11> + cz_sv = run_sv( + 'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[2] q;\nx q[0];\nx q[1];\ncz q[0], q[1];\n' + ) + assert_sv_close( + cz_sv, reference_statevector([("x", [0], []), ("x", [1], []), ("cz", [0, 1], [])], 2) + ) + # swap |q1=0,q0=1> (index 1) -> |q1=1,q0=0> (index 2) + swap_sv = run_sv( + 'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[2] q;\nx q[0];\nswap q[0], q[1];\n' + ) + assert np.isclose(swap_sv[2], 1.0) + # cx control=q0,target=q1 on |01> (index 1) -> |11> (index 3) + cx_sv = run_sv('OPENQASM 3;\ninclude "stdgates.inc";\nqubit[2] q;\nx q[0];\ncx q[0], q[1];\n') + assert np.isclose(cx_sv[3], 1.0) + + +def test_u3_fast_path_matches_decomposition(): + qasm = 'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[1] q;\nh q[0];\nu3(0.7, 1.1, -0.4) q[0];\n' + assert_sv_close(run_sv(qasm, external_gates=["u3"]), run_sv(qasm)) + + +# -------------------------------------------------------------------------- +# Gate fusion: a run of diagonal gates fuses into phases, and a following +# non-diagonal gate forces the accumulated matrix path + flush. +# -------------------------------------------------------------------------- + + +def test_diagonal_then_nondiagonal_fusion(): + qasm = ( + 'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[1] q;\n' + "h q[0];\n" # non-diagonal + "z q[0];\ns q[0];\nt q[0];\n" # fused diagonal run + "h q[0];\n" # forces flush + non-diagonal accumulation + "rz(0.5) q[0];\n" + ) + sv = run_sv(qasm) + expected = reference_statevector( + [ + ("h", [0], []), + ("z", [0], []), + ("s", [0], []), + ("t", [0], []), + ("h", [0], []), + ("rz", [0], [0.5]), + ], + 1, + ) + assert_sv_close(sv, expected) + + +def test_diagonal_gates_commute_through_cz(): + # Diagonal single-qubit gates remain fused across CZ; result must still be correct. + qasm = ( + 'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[2] q;\n' + "h q[0];\nh q[1];\nz q[0];\ncz q[0], q[1];\ns q[1];\n" + ) + sv = run_sv(qasm) + expected = reference_statevector( + [("h", [0], []), ("h", [1], []), ("z", [0], []), ("cz", [0, 1], []), ("s", [1], [])], + 2, + ) + assert_sv_close(sv, expected) + + +# -------------------------------------------------------------------------- +# Error handling and edge cases. +# -------------------------------------------------------------------------- + + +def test_module_not_unrolled_raises(): + """A QasmModule that was never unrolled must raise instead of silently + simulating the raw AST (which would drop control flow).""" + module = loads('OPENQASM 3;\ninclude "stdgates.inc";\nqubit[1] q;\nsx q[0];\n') + with pytest.raises(ValueError, match="has not been unrolled"): + Simulator().run(module, shots=0) + + +def test_unsupported_single_qubit_gate_raises(): + """The unsupported-gate branch is defense in depth: every stdgates gate + unrolls to the supported basis, so force an unknown name into the AST.""" + module = loads('OPENQASM 3;\ninclude "stdgates.inc";\nqubit[1] q;\nh q[0];\n') + module.unroll() + for statement in module._unrolled_ast.statements: + if type(statement).__name__ == "QuantumGate": + statement.name.name = "bogus" + with pytest.raises(ValueError, match="not supported by the statevector simulator"): + Simulator().run(module, shots=0) + + +def test_unsupported_two_qubit_gate_raises(): + module = loads('OPENQASM 3;\ninclude "stdgates.inc";\nqubit[2] q;\ncx q[0], q[1];\n') + module.unroll() + for statement in module._unrolled_ast.statements: + if type(statement).__name__ == "QuantumGate": + statement.name.name = "bogus2q" + with pytest.raises(ValueError, match="not supported by the statevector simulator"): + Simulator().run(module, shots=0) + + +def test_negative_shots_raises(): + with pytest.raises(ValueError, match="Shots must be"): + Simulator().run('OPENQASM 3;\ninclude "stdgates.inc";\nqubit[1] q;\nh q[0];\n', shots=-1) + + +def test_empty_circuit_returns_ground_state(): + result = Simulator(seed=0).run("OPENQASM 3;\nqubit[2] q;\n", shots=5) + assert np.isclose(result.final_statevector[0], 1.0) + assert sum(result.measurement_counts.values()) == 5 + + +def test_shots_sampling_is_seed_deterministic(): + qasm = 'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[2] q;\nh q[0];\ncx q[0], q[1];\n' + first = Simulator(seed=7).run(qasm, shots=128).measurement_counts + second = Simulator(seed=7).run(qasm, shots=128).measurement_counts + assert first == second + assert sum(first.values()) == 128 + # Bell state only populates |00> and |11>; each key is num_qubits long. + assert set(first) <= {"00", "11"} + + +def test_string_input_removes_idle_qubits(): + # q[1] is idle, so the reported circuit collapses to a single qubit. + result = Simulator(seed=0).run( + 'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[2] q;\nx q[0];\n', shots=4 + ) + assert len(result.final_statevector) == 2 + assert set(result.measurement_counts) == {"1"} + + +# -------------------------------------------------------------------------- +# Helper-level unit tests. +# -------------------------------------------------------------------------- + + +def test_rz_helper_matrix(): + mat = rz(np.pi / 2) + expected = np.array([[np.exp(-1j * np.pi / 4), 0], [0, np.exp(1j * np.pi / 4)]]) + assert np.allclose(mat, expected) + + +def test_try_eval_expression_constants_and_operators(): + assert _try_eval_expression(IntegerLiteral(value=3)) == 3.0 + assert _try_eval_expression(FloatLiteral(value=1.5)) == 1.5 + assert _try_eval_expression(Identifier(name="pi")) == pytest.approx(np.pi) + assert _try_eval_expression(Identifier(name="tau")) == pytest.approx(2 * np.pi) + assert _try_eval_expression(Identifier(name="euler")) == pytest.approx(np.e) + # unknown identifier is not evaluable + assert _try_eval_expression(Identifier(name="theta")) is None + + def binop(op, lhs, rhs): + return BinaryExpression( + op=BinaryOperator[op], lhs=FloatLiteral(value=lhs), rhs=FloatLiteral(value=rhs) + ) + + assert _try_eval_expression(binop("+", 1, 2)) == 3.0 + assert _try_eval_expression(binop("-", 5, 2)) == 3.0 + assert _try_eval_expression(binop("*", 3, 4)) == 12.0 + assert _try_eval_expression(binop("/", 8, 2)) == 4.0 + assert _try_eval_expression(binop("**", 2, 3)) == 8.0 + # unevaluable operand propagates as None + nested = BinaryExpression( + op=BinaryOperator["+"], lhs=Identifier(name="x"), rhs=FloatLiteral(value=1.0) + ) + assert _try_eval_expression(nested) is None + + neg = UnaryExpression(op=UnaryOperator["-"], expression=FloatLiteral(value=2.5)) + assert _try_eval_expression(neg) == -2.5 + assert ( + _try_eval_expression( + UnaryExpression(op=UnaryOperator["-"], expression=Identifier(name="x")) + ) + is None + ) + + # Operators / node types the evaluator does not handle return None rather + # than raising, so an un-evaluable parameter is simply skipped downstream. + assert _try_eval_expression(binop("%", 5, 2)) is None + assert ( + _try_eval_expression( + UnaryExpression(op=UnaryOperator["!"], expression=FloatLiteral(value=1.0)) + ) + is None + ) + assert _try_eval_expression(BooleanLiteral(value=True)) is None + + +def test_zero_angle_rotation_is_identity(): + """A lone zero-angle rotation produces a matrix tagged non-diagonal that is + actually the (diagonal) identity; the fusion flush must detect and apply it + as a no-op. ``rx(0)`` is kept alone on q[0] so it is not fused with another + non-diagonal gate first.""" + sv = run_sv('OPENQASM 3;\ninclude "stdgates.inc";\nqubit[2] q;\nh q[1];\nrx(0) q[0];\n') + expected = reference_statevector([("h", [1], [])], 2) + assert_sv_close(sv, expected) + + +def loads_unrolled(qasm): + module = loads(qasm) + module.unroll() + return module + + +# -------------------------------------------------------------------------- +# Multi-register programs: operands must resolve through (register, index), +# not the register-local index alone. Regression for the silent wrong-answer +# bug where `a[0]` and `b[0]` both mapped to global qubit 0. +# -------------------------------------------------------------------------- + + +def test_multi_register_bell_state(): + """Entangling across two registers must correlate the right qubits.""" + module = loads_unrolled( + 'OPENQASM 3;\ninclude "stdgates.inc";\n' + "qubit[2] a;\nqubit[2] b;\nh a[0];\ncx a[0], b[0];\n" + ) + result = Simulator(seed=0).run(module, shots=0) + # a[0] is global qubit 0, b[0] is global qubit 2 -> |0000> + |0101> + expected = np.zeros(16) + expected[0b0000] = 0.5 + expected[0b0101] = 0.5 + assert np.allclose(result.probabilities, expected, atol=1e-12) + + +def test_multi_register_distinct_qubits(): + """x a[0]; h b[0] must act on two different qubits, not fuse onto one.""" + result = Simulator(seed=0).run( + 'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[1] a;\nqubit[1] b;\nx a[0];\nh b[0];\n', + shots=0, + ) + assert np.allclose(result.probabilities, [0, 0.5, 0, 0.5], atol=1e-12) + + +def test_register_index_out_of_range_raises(): + module = loads_unrolled('OPENQASM 3;\ninclude "stdgates.inc";\nqubit[2] q;\nh q[0];\n') + for statement in module._unrolled_ast.statements: + if type(statement).__name__ == "QuantumGate": + statement.qubits[0].indices[0][0].value = 5 + with pytest.raises(ValueError, match="out of range"): + Simulator().run(module, shots=0) + + +# -------------------------------------------------------------------------- +# ccx / toffoli: previously silently dropped (no else branch for 3-qubit +# gates), corrupting every circuit containing one. +# -------------------------------------------------------------------------- + + +def _reference_ccx(num_qubits, controls, target): + dim = 2**num_qubits + mat = np.eye(dim, dtype=complex) + i0 = (1 << controls[0]) | (1 << controls[1]) + i1 = i0 | (1 << target) + mat[i0, i0] = 0 + mat[i1, i1] = 0 + mat[i0, i1] = 1 + mat[i1, i0] = 1 + return mat + + +def test_ccx_truth_table(): + """ccx flips the target iff both controls are set.""" + for prep, expected_index in [ + ("", 0b000), + ("x q[0];", 0b001), + ("x q[1];", 0b010), + ("x q[0];\nx q[1];", 0b111), # both controls -> target flipped + ("x q[0];\nx q[1];\nx q[2];", 0b011), # target un-flipped + ]: + result = Simulator(seed=0).run( + 'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[3] q;\n' + f"{prep}\nccx q[0], q[1], q[2];\n", + shots=0, + ) + expected = np.zeros(8) + expected[expected_index] = 1 + assert np.allclose(result.probabilities, expected, atol=1e-12), prep + + +def test_ccx_statevector_exact(): + """ccx on a non-trivial superposition matches the exact unitary, phases included.""" + sv = run_sv( + 'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[3] q;\n' + "ry(0.7) q[0];\nry(1.1) q[1];\nry(0.4) q[2];\nccx q[0], q[1], q[2];\n" + ) + + def _ry(theta): + c, s = np.cos(theta / 2), np.sin(theta / 2) + return np.array([[c, -s], [s, c]], dtype=complex) + + prep = np.zeros(8, dtype=complex) + prep[0] = 1.0 + for qubit, theta in [(0, 0.7), (1, 1.1), (2, 0.4)]: + prep = _embed(_ry(theta), qubit, 3) @ prep + expected = _reference_ccx(3, (0, 1), 2) @ prep + assert np.allclose(sv, expected, atol=1e-12) + + +def test_ccx_nonadjacent_qubits(): + """ccx with non-adjacent and permuted operand order.""" + sv = run_sv( + 'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[4] q;\n' + "x q[3];\nx q[1];\nccx q[3], q[1], q[0];\n" + ) + expected = np.zeros(16, dtype=complex) + expected[0b1011] = 1 # q3, q1 set -> q0 flipped + assert np.allclose(sv, expected, atol=1e-12) + + +def test_ccx_identical_qubits_raises(): + module = loads_unrolled( + 'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[3] q;\nccx q[0], q[1], q[2];\n' + ) + for statement in module._unrolled_ast.statements: + if type(statement).__name__ == "QuantumGate": + statement.qubits[1].indices[0][0].value = 0 # duplicate control + with pytest.raises(ValueError, match="distinct"): + Simulator().run(module, shots=0) + + +# -------------------------------------------------------------------------- +# Statements with real semantics must raise instead of being silently +# dropped; harmless statements must still be accepted. +# -------------------------------------------------------------------------- + + +def test_reset_raises(): + with pytest.raises(NotImplementedError, match="reset"): + Simulator().run( + 'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[1] q;\nx q[0];\nreset q[0];\n', + shots=0, + ) + + +def test_mid_circuit_measurement_raises(): + with pytest.raises(NotImplementedError, match="[Mm]id-circuit"): + Simulator().run( + 'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[2] q;\nbit[2] c;\n' + "h q[0];\nc[0] = measure q[0];\nx q[1];\n", + shots=0, + ) + + +def test_branching_raises(): + module = loads( + 'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[2] q;\nbit[2] c;\n' + "h q[0];\nc[0] = measure q[0];\nif (c[0] == 1) { x q[1]; }\n" + ) + module.unroll() + with pytest.raises(NotImplementedError): + Simulator().run(module, shots=0) + + +def test_terminal_measurement_accepted(): + """Terminal measurement is honored by the sampling step and must not raise.""" + result = Simulator(seed=0).run( + 'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[1] q;\nbit[1] c;\n' + "h q[0];\nc[0] = measure q[0];\n", + shots=100, + ) + assert sum(result.measurement_counts.values()) == 100 + + +def test_barrier_is_noop(): + result = Simulator(seed=0).run( + 'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[2] q;\nh q[0];\nbarrier q;\ncx q[0], q[1];\n', + shots=0, + ) + assert np.allclose(result.probabilities, [0.5, 0, 0, 0.5], atol=1e-12) + + +# -------------------------------------------------------------------------- +# p / sx support: main's c3x/c4x decompositions emit literal `p` gates, which +# previously failed with a misleading "call unroll() first" error. +# -------------------------------------------------------------------------- + + +def test_p_gate_exact(): + """h; p(theta); h on |0>: amplitudes ((1+e^it)/2, (1-e^it)/2). + + ``external_gates=["p"]`` keeps the literal `p` in the AST so this exercises + the simulator's own diagonal fast path (the unroller's rz-chain + decomposition would only agree up to global phase). + """ + theta = 0.37 + sv = run_sv( + f'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[1] q;\nh q[0];\np({theta}) q[0];\nh q[0];\n', + external_gates=["p"], + ) + expected = np.array([(1 + np.exp(1j * theta)) / 2, (1 - np.exp(1j * theta)) / 2]) + assert np.allclose(sv, expected, atol=1e-12) + # And the decomposed path agrees up to global phase. + decomposed = run_sv( + f'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[1] q;\nh q[0];\np({theta}) q[0];\nh q[0];\n' + ) + assert_sv_close(decomposed, expected) + + +def test_sx_squared_is_x(): + sv = run_sv('OPENQASM 3;\ninclude "stdgates.inc";\nqubit[1] q;\nsx q[0];\nsx q[0];\n') + assert np.allclose(np.abs(sv) ** 2, [0, 1], atol=1e-12) + + +def test_c3x_runs(): + """c3x decomposes through `p` gates on main; the simulator must accept it.""" + result = Simulator(seed=0).run( + 'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[4] q;\n' + "x q[0];\nx q[1];\nx q[2];\nc3x q[0], q[1], q[2], q[3];\n", + shots=0, + ) + expected = np.zeros(16) + expected[0b1111] = 1 + assert np.allclose(result.probabilities, expected, atol=1e-9) + + +def test_c4x_runs(): + result = Simulator(seed=0).run( + 'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[5] q;\n' + "x q[0];\nx q[1];\nx q[2];\nx q[3];\nc4x q[0], q[1], q[2], q[3], q[4];\n", + shots=0, + ) + expected = np.zeros(32) + expected[0b11111] = 1 + assert np.allclose(result.probabilities, expected, atol=1e-9) + + +# -------------------------------------------------------------------------- +# Result conventions and resource guards. +# -------------------------------------------------------------------------- + + +def test_counts_keys_match_probabilities_indexing(): + """Counts keys are the binary rendering of the probabilities index (qubit 0 + rightmost, matching qiskit), so probabilities[int(key, 2)] is the outcome's + probability. Regression for the reversed-keys bug.""" + result = Simulator(seed=0).run( + 'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[3] q;\nid q[0];\nid q[1];\nx q[2];\n', + shots=10, + ) + assert dict(result.measurement_counts) == {"100": 10} + assert result.probabilities[int("100", 2)] == pytest.approx(1.0) + + +def test_max_qubits_guard(monkeypatch): + """Over-ceiling programs raise instead of dying on the 2**n allocation.""" + monkeypatch.setenv("PYQASM_SIM_MAX_QUBITS", "8") + qasm = 'OPENQASM 3;\ninclude "stdgates.inc";\nqubit[9] q;\n' + "".join( + f"h q[{i}];\n" for i in range(9) + ) + with pytest.raises(ValueError, match="ceiling"): + Simulator().run(qasm, shots=0) + + +def test_shots_validated_before_simulation(): + with pytest.raises(ValueError, match="Shots must be"): + Simulator().run("OPENQASM 3;\nqubit[1] q;\n", shots=-1) + + +# -------------------------------------------------------------------------- +# Kernel argument validation: malformed calls raise ValueError instead of +# corrupting memory / crashing the interpreter. +# -------------------------------------------------------------------------- + + +def test_kernel_validation_rejects_malformed_calls(): + # pylint: disable-next=import-outside-toplevel,no-name-in-module + from pyqasm.accelerate import sv_sim + + identity = np.array([1, 0, 0, 1], dtype=np.complex128) + sv = np.zeros(8, dtype=np.complex128) + sv[0] = 1 + + with pytest.raises(ValueError, match="out of range"): + sv_sim.apply_single_qubit_gate(sv, 3, 40, identity) + with pytest.raises(ValueError, match="does not match"): + sv_sim.apply_single_qubit_gate(sv, 30, 0, identity) + with pytest.raises(ValueError, match="4 elements"): + sv_sim.apply_single_qubit_gate(sv, 3, 0, np.array([1], dtype=np.complex128)) + with pytest.raises(ValueError, match="both"): + sv_sim.apply_controlled_gate(sv, 3, 1, 1, identity) + with pytest.raises(ValueError, match="does not match"): + sv_sim.apply_single_qubit_gate(np.zeros(0, dtype=np.complex128), 1, 0, identity) + + +def test_apply_circuit_validation(): + # pylint: disable-next=import-outside-toplevel,no-name-in-module + from pyqasm.accelerate import sv_sim + + sv = np.zeros(8, dtype=np.complex128) + sv[0] = 1 + zeros_i32 = np.zeros(1, dtype=np.int32) + gate_params = np.zeros(4, dtype=np.complex128) + diag_phases = np.zeros(2, dtype=np.complex128) + tq_gates = np.zeros(16, dtype=np.complex128) + + with pytest.raises(ValueError, match="shorter than"): + sv_sim.apply_circuit( + sv, + 3, + zeros_i32, + zeros_i32, + zeros_i32, + gate_params, + diag_phases, + zeros_i32, + tq_gates, + 10_000_000, + ) + with pytest.raises(ValueError, match="out of range"): + sv_sim.apply_circuit( + sv, + 3, + np.array([0], dtype=np.int32), + np.array([40], dtype=np.int32), + zeros_i32, + gate_params, + diag_phases, + zeros_i32, + tq_gates, + 1, + ) + with pytest.raises(ValueError, match="offset"): + sv_sim.apply_circuit( + sv, + 3, + np.array([4], dtype=np.int32), + np.array([1], dtype=np.int32), + zeros_i32, + gate_params, + diag_phases, + np.array([1 << 28], dtype=np.int32), + tq_gates, + 1, + ) diff --git a/tox.ini b/tox.ini index 5aa7bca3..a00d624f 100644 --- a/tox.ini +++ b/tox.ini @@ -18,7 +18,7 @@ extras = cli pulse commands = - pytest tests --cov=pyqasm --cov-config=pyproject.toml --cov-report=term --cov-report=xml {posargs} + pytest tests -m "not benchmark" --cov=pyqasm --cov-config=pyproject.toml --cov-report=term --cov-report=xml {posargs} [testenv:docs] description = Use sphinx to build the HTML docs.