From ac30848e09598d9d9bba3f83c292badcce70ef15 Mon Sep 17 00:00:00 2001 From: Dmytro Fedoriaka Date: Sun, 27 Sep 2026 17:25:50 -0700 Subject: [PATCH 01/11] Add sparse simulator and tests --- .github/workflows/ci.yml | 6 +- cirq_sparse_sim/README.md | 33 + cirq_sparse_sim/__init__.py | 0 cirq_sparse_sim/sparse_sim.py | 436 +++++++++++ cirq_sparse_sim/sparse_sim_test.py | 1071 ++++++++++++++++++++++++++++ 5 files changed, 1543 insertions(+), 3 deletions(-) create mode 100644 cirq_sparse_sim/README.md create mode 100644 cirq_sparse_sim/__init__.py create mode 100644 cirq_sparse_sim/sparse_sim.py create mode 100644 cirq_sparse_sim/sparse_sim_test.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index c41af79..da63d3d 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -11,10 +11,10 @@ jobs: name: Tests (Ubuntu) runs-on: ubuntu-latest steps: - - uses: actions/checkout@v1 - - uses: actions/setup-python@v1 + - uses: actions/checkout@v7 + - uses: actions/setup-python@v7 with: - python-version: '3.12' + python-version: '3.14' architecture: 'x64' - name: Install requirements run: | diff --git a/cirq_sparse_sim/README.md b/cirq_sparse_sim/README.md new file mode 100644 index 0000000..bdb982b --- /dev/null +++ b/cirq_sparse_sim/README.md @@ -0,0 +1,33 @@ +# Sparse simulator for Cirq + +This is a simple implementation of a Sparse simulator for circuits in Cirq. + +It is sparse in a sense that it maintains a list of basis states with only +non-zero amplitudes. + +It supports unitary gates and computational basis measurments, and designed +primarily for simulating arithmetic circuits. + +I was not able to find such a simulator in Cirq (maybe I didn't look had +enough), so I just wrote this one for my needs. + +AI was used while writing this simulator. + +## Tests + +From the repository root, with Cirq, NumPy, SymPy, and pytest installed: + +```sh +python -m pytest cirq_sparse_sim/sparse_sim_test.py +``` + +The tests cover measurement results and final states, seeded random circuits +against Cirq's dense and Clifford simulators, arithmetic gates, classical +controls, reset, decomposition, qubit management, and error cases. Dense state +comparisons retain global phase; Clifford comparisons allow an overall global +phase because stabilizer representations need not preserve it. + +This suite targets the existing single-shot `run(circuit)` API and its +`cirq.Result` output. Full `cirq.SimulatorBase` compatibility, including +`simulate`, repetitions, parameter sweeps, and arbitrary qubit types, is not +implemented or claimed by these tests. \ No newline at end of file diff --git a/cirq_sparse_sim/__init__.py b/cirq_sparse_sim/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/cirq_sparse_sim/sparse_sim.py b/cirq_sparse_sim/sparse_sim.py new file mode 100644 index 0000000..484b8cb --- /dev/null +++ b/cirq_sparse_sim/sparse_sim.py @@ -0,0 +1,436 @@ +"""Sparse, single-shot simulator for a restricted subset of Cirq circuits. + +The simulator keeps the quantum state as a sparse list of computational basis +states and their amplitudes instead of a dense state vector. This is efficient +for circuits that mostly preserve sparsity (for example arithmetic/reversible +circuits with occasional single-qubit superposition gates). + +Execution model: + +1. Start in ``|0...0>`` represented as one basis state with amplitude ``1``. +2. Apply operations one by one, updating sparse basis states/amplitudes. +3. Measurements collapse the sparse state stochastically according to Born rule. +4. Return a single-shot ``cirq.ResultDict``. + +This simulator is intentionally limited and aimed at tests and debugging. +Unsupported operations raise ``ValueError``. +""" + +import math +import random +from collections.abc import Callable, Iterable, Sequence + +import cirq +import numpy as np +import sympy + +_AMPLITUDE_EPS = 1e-9 + + +class SparseSimArithmeticGate(cirq.Gate): + """Test-only gate that represents a reversible arithmetic basis-state map. + + States are interpreted as one or more little-endian integer registers. + ``f`` maps those integer register values to new register values. + ``f`` must represent a reversible transformation, but reversibility is not + checked by this class. + + Can only be simulated with SparseSimulator. + """ + + def __init__( + self, + input_sizes: list[int], + f: Callable[..., tuple[int, ...]], + ) -> None: + """Initialize the gate.""" + self.input_sizes = input_sizes + self.f = f + + def num_qubits(self) -> int: + """The number of qubits this gate acts on.""" + return sum(self.input_sizes) + + +class SparseSimQubitManager(cirq.QubitManager): + """Simple qubit manager with allocation/release over ``LineQubit`` indices.""" + + def __init__(self, simulator: "SparseSimulator") -> None: + """Initialize the simulator.""" + self.simulator = simulator + self.num_qubits = 0 + self.free_qubits: list[int] = [] + self.free_qubits_set: set[int] = set() + + def _allocate_qubit(self) -> cirq.LineQubit: + if len(self.free_qubits) == 0: + qubit_id = self.num_qubits + self.num_qubits += 1 + else: + qubit_id = self.free_qubits.pop() + self.free_qubits_set.remove(qubit_id) + return cirq.LineQubit(qubit_id) + + def qalloc(self, n: int, dim: int = 2) -> list[cirq.Qid]: + """Allocate qubits.""" + assert dim == 2 + return [self._allocate_qubit() for _ in range(n)] + + def qborrow(self, n: int, dim: int = 2) -> list[cirq.Qid]: + """Not implemented.""" + raise NotImplementedError("qborrow is not supported") + + def qfree(self, qubits: Iterable[cirq.Qid]) -> None: + """Free qubits.""" + for q in qubits: + assert isinstance(q, cirq.LineQubit) + assert 0 <= q.x < self.num_qubits + if q.x in self.free_qubits_set: + raise RuntimeError(f"Qubit {q.x} released twice") + if not self.simulator._is_qubit_zero(q.x): + raise RuntimeError(f"Qubit {q.x} released not in zero state") + self.free_qubits.append(q.x) + self.free_qubits_set.add(q.x) + + def num_allocated_qubits(self) -> int: + """Returns number of allocated qubits.""" + return self.num_qubits - len(self.free_qubits) + + +def _apply_cx(state: int, control_bit: int, target_bit: int) -> int: + if (state >> control_bit) % 2 == 1: + return state ^ (1 << target_bit) + return state + + +def _apply_ccx(state: int, ctrl1: int, ctrl2: int, target: int) -> int: + if (state >> ctrl1) % 2 == 1 and (state >> ctrl2) % 2 == 1: + return state ^ (1 << target) + return state + + +def _cz_phase(state: int, ctrl: int, target: int, exponent: float) -> complex: + if (state >> ctrl) % 2 == 1 and (state >> target) % 2 == 1: + return complex(math.cos(math.pi * exponent), math.sin(math.pi * exponent)) + return 1.0 + 0.0j + + +def _normalize(amplitudes: list[complex]) -> list[complex]: + norm = math.sqrt(sum(abs(a) ** 2 for a in amplitudes)) + return [a / norm for a in amplitudes] + + +class SparseSimulator: + """Sparse simulator for a restricted set of Cirq operations. + + Supported operations include: + - ``X``, ``CNOT``, ``CCNOT`` (basis permutation). + - ``CZPowGate`` with numeric exponent (phase on ``|11>`` branch). + - ``SparseSimArithmeticGate``. + - ``MeasurementGate`` and ``ResetChannel``. + - ``ClassicallyControlledOperation`` with ``KeyCondition`` or ``SympyCondition``. + - Any single-qubit gate with known unitary. + + Other operations are attempted via one-step decomposition and processed + recursively. If an operation still cannot be handled, simulation fails. + + Usage: + - Create a simulator instance: ``sim = SparseSimulator()``. + - Allocate qubits using ``sim.qubit_manager.qalloc``. Qubits created in any other + way are not supported. + - Create a circuit. + - Call ``sim.run(circuit)``. + - Use ``sim.read_register`` to read the value stored in a given register, + interpreted as an unsigned little-endian integer. + - Inspect measurement results in the ``cirq.Result`` object returned by + ``sim.run``. + + Notes: + - ``run`` is single-shot only. + - ``measurement_results`` stores the latest integer value per measurement key. + These values are little-endian; classical conditions use Cirq's big-endian + integers and measurement occurrence indices. + - ``read_register`` requires the final state to be a single basis state. + If it's not the case, add measurements. + + """ + + def __init__(self, seed: int | None = None) -> None: + """Initialize SparseSimulator.""" + self.random = random.Random(seed) + self.qubit_manager = SparseSimQubitManager(self) + + # Sparse state is represented by a list of basis states with + # non-zero amplitudes and a list of amplitudes. + self.basis_states = [0] + self.amplitudes = [1.0 + 0.0j] + + # Measurement results. + self.measurement_results: dict[str, int] = {} + self._measurement_history: dict[str, list[int]] = {} + + # Length of measurements (to restore measurements as bit vectors). + self._meas_len: dict[str, int] = {} + + def _apply_x(self, qid: int) -> None: + self.basis_states = [s ^ (1 << qid) for s in self.basis_states] + + def _apply_cx(self, ctrl: int, target: int) -> None: + self.basis_states = [_apply_cx(s, ctrl, target) for s in self.basis_states] + + def _apply_ccx(self, ctrl1: int, ctrl2: int, target: int) -> None: + self.basis_states = [ + _apply_ccx(s, ctrl1, ctrl2, target) for s in self.basis_states + ] + + def _apply_cz(self, ctrl: int, target: int, exponent: float) -> None: + self.amplitudes = [ + amplitude * _cz_phase(state, ctrl, target, exponent) + for state, amplitude in zip(self.basis_states, self.amplitudes, strict=True) + ] + + def _apply_global_shift(self, gate: cirq.EigenGate) -> None: + if gate.global_shift == 0: + return + angle = math.pi * float(gate.exponent) * gate.global_shift + phase = complex(math.cos(angle), math.sin(angle)) + self.amplitudes = [amplitude * phase for amplitude in self.amplitudes] + + def _apply_arithmetic_gate( + self, op_gate: SparseSimArithmeticGate, qubits: list[int] + ) -> None: + assert len(qubits) == op_gate.num_qubits() + registers: list[list[int]] = [] + offset = 0 + for n in op_gate.input_sizes: + registers.append(qubits[offset : offset + n]) + offset += n + new_basis_states = [] + for s in self.basis_states: + inputs = [ + sum(((s >> q) & 1) << i for i, q in enumerate(r)) for r in registers + ] + outputs = op_gate.f(*inputs) + assert len(outputs) == len(inputs) + new_state = s + for inp, out, reg in zip(inputs, outputs, registers, strict=True): + xor_val = inp ^ out + for i, q in enumerate(reg): + if (xor_val >> i) & 1: + new_state ^= 1 << q + new_basis_states.append(new_state) + self.basis_states = new_basis_states + + def _apply_single_qubit_unitary_gate( + self, op_gate: cirq.Gate, qubit_id: int + ) -> None: + """Apply a 2x2 unitary to one qubit across all sparse basis states. + + Output contributions that land on the same basis state are summed. + Terms with very small amplitudes (``<_AMPLITUDE_EPS``) are dropped. + """ + matrix = cirq.unitary(op_gate, default=None) + assert matrix is not None + assert matrix.shape == (2, 2) + + merged_amplitudes: dict[int, complex] = {} + for state, amplitude in zip(self.basis_states, self.amplitudes, strict=True): + input_bit = (state >> qubit_id) & 1 + for output_bit in (0, 1): + coeff = matrix[output_bit, input_bit] + if coeff == 0: + continue + new_state = ( + state if output_bit == input_bit else state ^ (1 << qubit_id) + ) + merged_amplitudes[new_state] = ( + merged_amplitudes.get(new_state, 0.0 + 0.0j) + amplitude * coeff + ) + + kept_states = [ + (state, amplitude) + for state, amplitude in merged_amplitudes.items() + if abs(amplitude) >= _AMPLITUDE_EPS + ] + assert len(kept_states) > 0 + self.basis_states = [state for state, _ in kept_states] + self.amplitudes = _normalize([amplitude for _, amplitude in kept_states]) + + def _measure_qubit(self, target_qubit: int) -> int: + """Measure one qubit and collapse sparse state according to Born rule. + + Returns the measured bit (0/1). If all current basis states already + agree on the measured bit, this is a cheap no-op collapse. + """ + if len(self.basis_states) == 1: + return (self.basis_states[0] >> target_qubit) & 1 + + zero_subspace: list[tuple[int, complex]] = [] + one_subspace: list[tuple[int, complex]] = [] + for state, amplitude in zip(self.basis_states, self.amplitudes, strict=True): + if (state >> target_qubit) & 1: + one_subspace.append((state, amplitude)) + else: + zero_subspace.append((state, amplitude)) + + if len(one_subspace) == 0: + return 0 + if len(zero_subspace) == 0: + return 1 + + zero_prob = sum(abs(amplitude) ** 2 for _, amplitude in zero_subspace) + one_prob = sum(abs(amplitude) ** 2 for _, amplitude in one_subspace) + assert zero_prob > 0 + assert one_prob > 0 + + if self.random.random() < zero_prob / (zero_prob + one_prob): + measured_bit = 0 + chosen = zero_subspace + else: + measured_bit = 1 + chosen = one_subspace + + self.basis_states = [state for state, _ in chosen] + self.amplitudes = _normalize([amplitude for _, amplitude in chosen]) + return measured_bit + + def _apply_measurement_gate( + self, op: cirq.Operation, op_gate: cirq.MeasurementGate, qubit_ids: list[int] + ) -> None: + """Apply a measurement gate and record key result as little-endian int.""" + key = cirq.measurement_key_name(op) + measured_value = 0 + measured_bits = [] + for i, qid in enumerate(qubit_ids): + measured_bit = self._measure_qubit(qid) + invert = i < len(op_gate.invert_mask) and op_gate.invert_mask[i] + if invert: + measured_bit ^= 1 + measured_value |= measured_bit << i + measured_bits.append(measured_bit) + self.measurement_results[key] = measured_value + self._meas_len[key] = len(qubit_ids) + self._measurement_history.setdefault(key, []).append( + cirq.big_endian_bits_to_int(measured_bits) + ) + + def _apply_reset(self, target_qubit: int) -> None: + """Reset one qubit to ``|0>`` via measure-and-conditional-flip.""" + if self._measure_qubit(target_qubit) == 1: + self._apply_x(target_qubit) + + def _resolve_sympy_condition(self, condition: cirq.SympyCondition) -> bool: + replacements: dict[str, object] = {} + for symbol in condition.expr.free_symbols: + if isinstance(symbol, sympy.Symbol): + replacements[symbol.name] = self._measurement_history[symbol.name][-1] + value = condition.expr.subs(replacements) + assert len(value.free_symbols) == 0 + return bool(value) + + def _classical_condition_satisfied(self, condition: cirq.Condition) -> bool: + """Evaluate a supported classical condition using stored measurements.""" + if isinstance(condition, cirq.SympyCondition): + return self._resolve_sympy_condition(condition) + if isinstance(condition, cirq.KeyCondition): + return self._measurement_history[str(condition.key)][condition.index] != 0 + raise ValueError(f"Unsupported classical condition: {condition}") + + def _run_op(self, op: cirq.Operation, context: cirq.DecompositionContext) -> None: + """Execute one operation, recursively decomposing when needed.""" + if isinstance(op, cirq.ClassicallyControlledOperation): + if all( + self._classical_condition_satisfied(condition) + for condition in op.classical_controls + ): + self._run_op(op.without_classical_controls(), context) + return + + op_gate = op.gate + qubit_ids = [q.x for q in op.qubits] + + if isinstance(op_gate, cirq.XPowGate) and op_gate.exponent == 1: + self._apply_x(qubit_ids[0]) + self._apply_global_shift(op_gate) + return + if isinstance(op_gate, cirq.CXPowGate) and op_gate.exponent == 1: + self._apply_cx(qubit_ids[0], qubit_ids[1]) + self._apply_global_shift(op_gate) + return + if isinstance(op_gate, cirq.CCXPowGate) and op_gate.exponent == 1: + self._apply_ccx(qubit_ids[0], qubit_ids[1], qubit_ids[2]) + self._apply_global_shift(op_gate) + return + if isinstance(op_gate, cirq.CZPowGate) and not cirq.is_parameterized(op_gate): + self._apply_cz(qubit_ids[0], qubit_ids[1], float(op_gate.exponent)) + self._apply_global_shift(op_gate) + return + if isinstance(op_gate, SparseSimArithmeticGate): + self._apply_arithmetic_gate(op_gate, qubit_ids) + return + if isinstance(op_gate, cirq.MeasurementGate): + self._apply_measurement_gate(op, op_gate, qubit_ids) + return + if isinstance(op_gate, cirq.ResetChannel): + assert len(qubit_ids) == 1 + self._apply_reset(qubit_ids[0]) + return + + # Prefer explicit unitaries to decompositions introducing unsupported gates. + if ( + op_gate is not None + and len(qubit_ids) == 1 + and cirq.has_unitary(op_gate, allow_decompose=False) + ): + self._apply_single_qubit_unitary_gate(op_gate, qubit_ids[0]) + return + + num_alloc_before = self.qubit_manager.num_allocated_qubits() + decomposed = cirq.decompose_once( + op, default=None, context=context, flatten=False + ) + if decomposed is not None: + for sub_op in cirq.flatten_to_ops(decomposed): + self._run_op(sub_op, context) + if self.qubit_manager.num_allocated_qubits() != num_alloc_before: + raise RuntimeError(f"Operation {op} did not free all allocated qubits") + return + + raise ValueError(f"Operation cannot be simulated: {op}.") + + def run(self, circuit: cirq.Circuit) -> cirq.Result: + """Run ``circuit`` once and return ``cirq.ResultDict``. + + Returned ``measurements[key]`` is a boolean array with shape + ``(1, num_measured_qubits_for_key)``. + """ + self.basis_states = [0] + self.amplitudes = [1.0 + 0.0j] + self.measurement_results = {} + self._measurement_history = {} + self._meas_len = {} + # Circuit must use only qubits allocated using qubit manager of this simulator. + context = cirq.DecompositionContext(self.qubit_manager) + for op in circuit.all_operations(): + self._run_op(op, context) + + measurements = { + key: np.asarray( + [[((value >> i) & 1) == 1 for i in range(self._meas_len[key])]], + dtype=np.bool_, + ) + for key, value in self.measurement_results.items() + } + return cirq.ResultDict(params=cirq.ParamResolver({}), measurements=measurements) + + def read_register(self, register: Sequence[cirq.Qid]) -> int: + """Read a little-endian register as an integer.""" + assert len(self.basis_states) == 1, ( + "Final state is superposition, add measurements." + ) + s = self.basis_states[0] + return sum(((s >> q.x) & 1) << i for i, q in enumerate(register)) + + def _is_qubit_zero(self, qid: int) -> bool: + """Checks whether qubit is in 0 state.""" + return all((s >> qid) & 1 == 0 for s in self.basis_states) \ No newline at end of file diff --git a/cirq_sparse_sim/sparse_sim_test.py b/cirq_sparse_sim/sparse_sim_test.py new file mode 100644 index 0000000..3ad5f30 --- /dev/null +++ b/cirq_sparse_sim/sparse_sim_test.py @@ -0,0 +1,1071 @@ +"""Tests for SparseSimulator's existing single-shot API.""" + +import itertools +import math +from collections.abc import Iterator, Sequence + +import cirq +import numpy as np +import pytest +import sympy + +from cirq_sparse_sim.sparse_sim import SparseSimArithmeticGate, SparseSimulator + + +def _state_vector( + simulator: SparseSimulator, qubit_order: Sequence[cirq.Qid] +) -> np.ndarray: + """Convert little-endian sparse indices to Cirq's specified tensor order.""" + qubit_ids = [] + for qubit in qubit_order: + assert isinstance(qubit, cirq.LineQubit) + qubit_ids.append(qubit.x) + + assert len(simulator.basis_states) == len(simulator.amplitudes) + assert len(set(simulator.basis_states)) == len(simulator.basis_states) + assert np.all(np.isfinite(simulator.amplitudes)) + np.testing.assert_allclose( + np.linalg.norm(simulator.amplitudes), 1, rtol=0, atol=1e-12 + ) + + vector = np.zeros(1 << len(qubit_ids), dtype=np.complex128) + included_bits = sum(1 << qid for qid in qubit_ids) + for basis, amplitude in zip( + simulator.basis_states, simulator.amplitudes, strict=True + ): + assert basis >= 0 + assert basis & ~included_bits == 0, "An omitted qubit is not in |0>." + index = 0 + for qid in qubit_ids: + index = (index << 1) | ((basis >> qid) & 1) + vector[index] = amplitude + return vector + + +def _prepare_basis( + qubits: Sequence[cirq.Qid], value: int +) -> list[cirq.Operation]: + return [cirq.X(q) for i, q in enumerate(qubits) if (value >> i) & 1] + + +def _assert_matches_dense( + simulator: SparseSimulator, + circuit: cirq.Circuit, + qubit_order: Sequence[cirq.Qid], +) -> cirq.Result: + """Compare circuits with no measurements or only deterministic outcomes.""" + actual = simulator.run(circuit) + expected = cirq.Simulator(dtype=np.complex128).simulate( + circuit, qubit_order=qubit_order + ) + np.testing.assert_allclose( + _state_vector(simulator, qubit_order), + expected.final_state_vector, + rtol=0, + atol=1e-8, + ) + assert actual.measurements.keys() == expected.measurements.keys() + for key, bits in expected.measurements.items(): + np.testing.assert_array_equal(actual.measurements[key], bits[np.newaxis, :]) + return actual + + +@pytest.mark.parametrize("num_qubits", [0, 1, 4]) +def test_empty_circuit(num_qubits: int) -> None: + simulator = SparseSimulator(seed=1) + qubits = simulator.qubit_manager.qalloc(num_qubits) + + result = _assert_matches_dense(simulator, cirq.Circuit(), qubits) + + assert isinstance(result, cirq.ResultDict) + assert result.params == cirq.ParamResolver({}) + assert result.measurements == {} + assert simulator.measurement_results == {} + assert simulator.basis_states == [0] + assert simulator.amplitudes == [1] + assert simulator.read_register(qubits) == 0 + assert simulator.read_register([]) == 0 + + +@pytest.mark.parametrize("initial", range(8)) +@pytest.mark.parametrize("gate", [cirq.X, cirq.CNOT, cirq.CCNOT], ids=repr) +def test_basis_permutations_measurements_and_final_state( + initial: int, gate: cirq.Gate +) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(3) + if gate == cirq.X: + operation = gate(qubits[2]) + expected = initial ^ 4 + elif gate == cirq.CNOT: + operation = gate(qubits[0], qubits[2]) + expected = initial ^ (4 if initial & 1 else 0) + else: + operation = gate(qubits[1], qubits[0], qubits[2]) + expected = initial ^ (4 if initial & 3 == 3 else 0) + measurement_order = [qubits[2], qubits[0], qubits[1]] + circuit = cirq.Circuit( + _prepare_basis(qubits, initial), + operation, + cirq.measure(*measurement_order, key="result"), + ) + + result = _assert_matches_dense(simulator, circuit, qubits) + + expected_bits = [(expected >> i) & 1 for i in (2, 0, 1)] + np.testing.assert_array_equal(result.measurements["result"], [expected_bits]) + assert simulator.basis_states == [expected] + assert simulator.read_register(qubits) == expected + assert simulator.read_register(qubits[::-1]) == cirq.big_endian_bits_to_int( + [(expected >> i) & 1 for i in range(3)] + ) + + +@pytest.mark.parametrize("initial", [0, 1]) +@pytest.mark.parametrize( + "gate", + [ + cirq.I, + cirq.X, + cirq.Y, + cirq.Z, + cirq.H, + cirq.S, + cirq.T, + cirq.S**-1, + cirq.T**-1, + cirq.X**0.25, + cirq.Y**-0.5, + cirq.Z**1.25, + cirq.rx(0.731), + cirq.ry(-0.413), + cirq.rz(1.271), + cirq.rx(np.pi), + cirq.XPowGate(exponent=1, global_shift=0.25), + pytest.param( + cirq.MatrixGate(cirq.testing.random_unitary(2, random_state=123)), + id="random-unitary", + ), + ], + ids=repr, +) +def test_single_qubit_unitaries(initial: int, gate: cirq.Gate) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(1) + circuit = cirq.Circuit(_prepare_basis(qubits, initial), gate(*qubits)) + + simulator.run(circuit) + + np.testing.assert_allclose( + _state_vector(simulator, qubits), + cirq.unitary(gate)[:, initial], + rtol=0, + atol=1e-9, + ) + + +@pytest.mark.parametrize( + "gate", + [ + cirq.CZ**0, + cirq.CZ**0.25, + cirq.CZ**-0.75, + cirq.CZ**2, + cirq.CZPowGate(exponent=0.5, global_shift=-0.5), + cirq.CXPowGate(exponent=1, global_shift=-0.5), + cirq.CCXPowGate(exponent=1, global_shift=0.25), + cirq.CNOT**0.5, + cirq.CCNOT**0.5, + cirq.SWAP, + cirq.ISWAP**0.5, + cirq.CCZ**0.25, + cirq.FSimGate(theta=0.3, phi=-0.7), + ], + ids=repr, +) +def test_multiqubit_gates_and_decompositions(gate: cirq.Gate) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(cirq.num_qubits(gate)) + circuit = cirq.Circuit( + [cirq.ry(0.2 + 0.3 * i)(q) for i, q in enumerate(qubits)], + cirq.T(qubits[0]), + gate(*qubits[::-1]), + ) + + _assert_matches_dense(simulator, circuit, qubits) + + +@pytest.mark.parametrize("initial", [0, 1]) +@pytest.mark.parametrize( + "gate", cirq.SingleQubitCliffordGate.all_single_qubit_cliffords, ids=repr +) +def test_all_single_qubit_cliffords(initial: int, gate: cirq.Gate) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(1) + + _assert_matches_dense( + simulator, + cirq.Circuit(_prepare_basis(qubits, initial), gate(*qubits)), + qubits, + ) + + +@pytest.mark.parametrize("seed", range(10)) +@pytest.mark.parametrize("num_qubits", [1, 2, 3, 5]) +def test_random_circuits_match_dense_simulator(seed: int, num_qubits: int) -> None: + simulator = SparseSimulator(seed=seed) + qubits = simulator.qubit_manager.qalloc(num_qubits) + gate_domain = { + cirq.X: 1, + cirq.Y: 1, + cirq.Z: 1, + cirq.H: 1, + cirq.S: 1, + cirq.T: 1, + cirq.X**0.37: 1, + cirq.Y**-0.23: 1, + cirq.rz(0.71): 1, + cirq.MatrixGate(cirq.testing.random_unitary(2, random_state=seed)): 1, + cirq.CNOT: 2, + cirq.CZ**0.31: 2, + cirq.SWAP: 2, + cirq.ISWAP**0.5: 2, + cirq.CCNOT: 3, + } + circuit = cirq.testing.random_circuit( + qubits, + n_moments=30, + op_density=0.85, + gate_domain=gate_domain, + random_state=seed, + ) + qubit_order = list(np.random.RandomState(seed).permutation(qubits)) + + _assert_matches_dense(simulator, circuit, qubit_order) + + +@pytest.mark.parametrize("seed", range(10)) +@pytest.mark.parametrize("num_qubits", [1, 2, 4, 6]) +def test_random_clifford_circuits_match_clifford_simulator( + seed: int, num_qubits: int +) -> None: + simulator = SparseSimulator(seed=seed) + qubits = simulator.qubit_manager.qalloc(num_qubits) + gate_domain = { + gate: 1 for gate in cirq.SingleQubitCliffordGate.all_single_qubit_cliffords + } + gate_domain.update({cirq.CNOT: 2, cirq.CZ: 2, cirq.SWAP: 2, cirq.ISWAP: 2}) + circuit = cirq.testing.random_circuit( + qubits, + n_moments=40, + op_density=0.85, + gate_domain=gate_domain, + random_state=seed, + ) + assert all(cirq.has_stabilizer_effect(op) for op in circuit.all_operations()) + qubit_order = qubits[::-1] + + simulator.run(circuit) + expected = cirq.CliffordSimulator(seed=seed).simulate( + circuit, qubit_order=qubit_order + ) + + # Stabilizer representations need not retain the circuit's global phase. + cirq.testing.assert_allclose_up_to_global_phase( + _state_vector(simulator, qubit_order), + expected.final_state.state_vector(), + atol=1e-7, + ) + + +@pytest.mark.parametrize("seed", range(5)) +def test_random_circuit_followed_by_inverse_returns_to_zero(seed: int) -> None: + simulator = SparseSimulator(seed=seed) + qubits = simulator.qubit_manager.qalloc(4) + circuit = cirq.testing.random_circuit( + qubits, n_moments=25, op_density=0.9, random_state=seed + ) + circuit += cirq.inverse(circuit) + circuit += cirq.measure(*qubits, key="zero") + + result = _assert_matches_dense(simulator, circuit, qubits) + + np.testing.assert_array_equal(result.measurements["zero"], [[0, 0, 0, 0]]) + assert simulator.basis_states == [0] + np.testing.assert_allclose(simulator.amplitudes, [1], atol=1e-9) + + +def test_destructive_interference_merges_and_prunes_states() -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(3) + circuit = cirq.Circuit( + cirq.H.on_each(*qubits), + cirq.Z(qubits[1]), + cirq.H.on_each(*qubits), + ) + + _assert_matches_dense(simulator, circuit, qubits) + + assert simulator.basis_states == [2] + np.testing.assert_allclose(simulator.amplitudes, [1], atol=1e-12) + + +def test_negligible_amplitudes_are_pruned_and_state_is_normalized() -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(1) + + simulator.run(cirq.Circuit(cirq.ry(1e-10)(qubits[0]))) + + assert simulator.basis_states == [0] + np.testing.assert_allclose(_state_vector(simulator, qubits), [1, 0], atol=1e-12) + + +def test_small_but_significant_amplitudes_are_retained() -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(1) + gate = cirq.ry(1e-8) + + simulator.run(cirq.Circuit(gate(*qubits))) + + assert set(simulator.basis_states) == {0, 1} + np.testing.assert_allclose( + _state_vector(simulator, qubits), cirq.unitary(gate)[:, 0], rtol=0, atol=1e-12 + ) + + +@pytest.mark.parametrize("seed", [0, 1, 2, 3]) +@pytest.mark.parametrize("num_qubits", [2, 3]) +def test_entangled_measurements_and_collapsed_phase( + seed: int, num_qubits: int +) -> None: + simulator = SparseSimulator(seed=seed) + qubits = simulator.qubit_manager.qalloc(num_qubits) + circuit = cirq.Circuit( + cirq.H(qubits[0]), + [cirq.CNOT(qubits[0], q) for q in qubits[1:]], + cirq.S(qubits[0]), + cirq.measure(qubits[-1], key="first"), + cirq.measure(*qubits, key="all"), + cirq.measure(qubits[0], key="again"), + ) + + result = simulator.run(circuit) + + bit = int(result.measurements["first"][0, 0]) + assert bit in (0, 1) + np.testing.assert_array_equal(result.measurements["all"], [[bit] * num_qubits]) + np.testing.assert_array_equal(result.measurements["again"], [[bit]]) + assert simulator.read_register(qubits) == bit * ((1 << num_qubits) - 1) + expected = np.zeros(1 << num_qubits, dtype=np.complex128) + expected[-1 if bit else 0] = 1j if bit else 1 + np.testing.assert_allclose(_state_vector(simulator, qubits), expected, atol=1e-12) + + +@pytest.mark.parametrize("seed", [0, 1]) +def test_partial_measurement_preserves_unmeasured_superposition(seed: int) -> None: + simulator = SparseSimulator(seed=seed) + qubits = simulator.qubit_manager.qalloc(3) + result = simulator.run( + cirq.Circuit( + cirq.H(qubits[0]), + cirq.CNOT(qubits[0], qubits[1]), + cirq.H(qubits[2]), + cirq.measure(qubits[0], key="m"), + ) + ) + + bit = int(result.measurements["m"][0, 0]) + expected = np.zeros(8, dtype=np.complex128) + expected[6 * bit : 6 * bit + 2] = 1 / math.sqrt(2) + np.testing.assert_allclose(_state_vector(simulator, qubits), expected, atol=1e-12) + with pytest.raises(AssertionError, match="superposition"): + simulator.read_register(qubits) + + +@pytest.mark.parametrize("bit", [0, 1]) +def test_deterministic_measurement_preserves_other_qubits(bit: int) -> None: + simulator = SparseSimulator(seed=0) + qubits = simulator.qubit_manager.qalloc(2) + circuit = cirq.Circuit( + _prepare_basis(qubits[:1], bit), + cirq.H(qubits[1]), + cirq.measure(qubits[0], key="m"), + ) + + result = _assert_matches_dense(simulator, circuit, qubits) + + np.testing.assert_array_equal(result.measurements["m"], [[bit]]) + assert len(simulator.basis_states) == 2 + + +@pytest.mark.parametrize("initial", [0, 5, 7]) +@pytest.mark.parametrize("invert_mask", [(), (True,), (True, False, True)]) +def test_measurement_invert_mask_changes_result_not_state( + initial: int, invert_mask: tuple[bool, ...] +) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(3) + measurement_order = [qubits[2], qubits[0], qubits[1]] + result = _assert_matches_dense( + simulator, + cirq.Circuit( + _prepare_basis(qubits, initial), + cirq.measure(*measurement_order, key="m", invert_mask=invert_mask), + ), + qubits, + ) + + bits = result.measurements["m"][0] + assert simulator.basis_states == [initial] + assert simulator.read_register(qubits) == initial + assert simulator.measurement_results["m"] == sum( + int(bit) << i for i, bit in enumerate(bits) + ) + + +@pytest.mark.parametrize("key", [None, "bits", cirq.MeasurementKey("bits", ("scope",))]) +def test_measurement_keys_and_result_interface( + key: str | cirq.MeasurementKey | None, +) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(3) + measurement = cirq.measure(*qubits, key=key) + name = cirq.measurement_key_name(measurement) + result = simulator.run( + cirq.Circuit( + cirq.X(qubits[0]), measurement, cirq.measure(qubits[0], key="flag") + ) + ) + + assert isinstance(result, cirq.Result) + assert result.params == cirq.ParamResolver({}) + assert result.measurements[name].shape == (1, 3) + assert result.measurements[name].dtype == np.bool_ + np.testing.assert_array_equal(result.measurements[name], [[1, 0, 0]]) + np.testing.assert_array_equal(result.records[name], [[[1, 0, 0]]]) + assert result.histogram(key=name) == {4: 1} + assert result.multi_measurement_histogram(keys=[name, "flag"]) == {(4, 1): 1} + assert result.data[name].tolist() == [4] + assert simulator.measurement_results[name] == 1 + assert simulator.read_register(qubits) == 1 + np.testing.assert_array_equal(_state_vector(simulator, qubits), np.eye(8)[4]) + + +def test_repeated_measurement_key_retains_latest_value() -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(2) + result = simulator.run( + cirq.Circuit( + cirq.X(qubits[0]), + cirq.measure(qubits[0], key="m"), + cirq.X(qubits[0]), + cirq.X(qubits[1]), + cirq.measure(*qubits, key="m"), + ) + ) + + assert simulator.measurement_results == {"m": 2} + np.testing.assert_array_equal(result.measurements["m"], [[0, 1]]) + assert simulator.read_register(qubits) == 2 + + +def test_measurements_follow_nonuniform_born_probabilities() -> None: + simulator = SparseSimulator(seed=12345) + qubits = simulator.qubit_manager.qalloc(1) + probability_one = 0.2 + trials = 1000 + circuit = cirq.Circuit( + cirq.ry(2 * math.asin(math.sqrt(probability_one)))(qubits[0]), + cirq.measure(qubits[0], key="m"), + ) + ones = 0 + + for _ in range(trials): + result = simulator.run(circuit) + bit = int(result.measurements["m"][0, 0]) + ones += bit + assert simulator.read_register(qubits) == bit + np.testing.assert_allclose(_state_vector(simulator, qubits), np.eye(2)[bit]) + + standard_deviation = math.sqrt(trials * probability_one * (1 - probability_one)) + assert abs(ones - trials * probability_one) < 6 * standard_deviation + + +def test_seed_reproducibility_across_repeated_runs() -> None: + outcomes = [] + for _ in range(2): + simulator = SparseSimulator(seed=2468) + qubits = simulator.qubit_manager.qalloc(2) + circuit = cirq.Circuit( + cirq.H.on_each(*qubits), cirq.measure(*qubits, key="m") + ) + outcomes.append( + [simulator.run(circuit).measurements["m"].tolist() for _ in range(32)] + ) + + assert outcomes[0] == outcomes[1] + assert len({tuple(shot[0]) for shot in outcomes[0]}) > 1 + + +def test_run_clears_previous_state_and_measurements_without_mutating_results() -> None: + simulator = SparseSimulator(seed=0) + qubits = simulator.qubit_manager.qalloc(2) + first = simulator.run( + cirq.Circuit(cirq.X.on_each(*qubits), cirq.measure(*qubits, key="old")) + ) + + second = simulator.run(cirq.Circuit(cirq.measure(qubits[0], key="new"))) + + np.testing.assert_array_equal(first.measurements["old"], [[1, 1]]) + np.testing.assert_array_equal(second.measurements["new"], [[0]]) + assert set(second.measurements) == {"new"} + assert simulator.measurement_results == {"new": 0} + assert simulator.basis_states == [0] + assert simulator.qubit_manager.num_allocated_qubits() == 2 + assert simulator.run(cirq.Circuit()).measurements == {} + assert simulator.measurement_results == {} + + +@pytest.mark.parametrize("initial", [0, 1]) +def test_reset_basis_state(initial: int) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(1) + result = _assert_matches_dense( + simulator, + cirq.Circuit( + _prepare_basis(qubits, initial), + cirq.reset(qubits[0]), + cirq.measure(qubits[0], key="m"), + ), + qubits, + ) + + np.testing.assert_array_equal(result.measurements["m"], [[0]]) + assert simulator.basis_states == [0] + + +@pytest.mark.parametrize("seed", [0, 1]) +def test_reset_entangled_qubit_collapses_partner_without_recording_result( + seed: int, +) -> None: + simulator = SparseSimulator(seed=seed) + qubits = simulator.qubit_manager.qalloc(2) + result = simulator.run( + cirq.Circuit( + cirq.H(qubits[0]), + cirq.CNOT(*qubits), + cirq.reset(qubits[0]), + cirq.measure(qubits[1], key="partner"), + ) + ) + + bit = int(result.measurements["partner"][0, 0]) + assert set(result.measurements) == {"partner"} + assert simulator.measurement_results == {"partner": bit} + assert simulator.read_register(qubits) == 2 * bit + np.testing.assert_allclose(_state_vector(simulator, qubits), np.eye(4)[bit]) + simulator.qubit_manager.qfree(qubits[:1]) + + +@pytest.mark.parametrize("bit", [0, 1]) +@pytest.mark.parametrize("invert", [False, True]) +def test_key_classical_control_uses_reported_measurement( + bit: int, invert: bool +) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(2) + circuit = cirq.Circuit( + _prepare_basis(qubits[:1], bit), + cirq.measure(qubits[0], key="control", invert_mask=(invert,)), + cirq.X(qubits[1]).with_classical_controls("control"), + cirq.measure(qubits[1], key="target"), + ) + + result = _assert_matches_dense(simulator, circuit, qubits) + + np.testing.assert_array_equal(result.measurements["target"], [[bit ^ invert]]) + assert simulator.read_register(qubits) == bit + 2 * (bit ^ invert) + + +@pytest.mark.parametrize("index", [-1, 0]) +def test_key_condition_checks_entire_register_not_indexed_bit(index: int) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(3) + condition = cirq.KeyCondition(cirq.MeasurementKey("control"), index=index) + circuit = cirq.Circuit( + cirq.X(qubits[1]), + cirq.measure(*qubits[:2], key="control"), + cirq.X(qubits[2]).with_classical_controls(condition), + cirq.measure(qubits[2], key="target"), + ) + + result = _assert_matches_dense(simulator, circuit, qubits) + + np.testing.assert_array_equal(result.measurements["target"], [[1]]) + + +@pytest.mark.parametrize("index, expected", [(0, 0), (1, 1), (-1, 1), (-2, 0)]) +def test_key_condition_index_selects_measurement_occurrence( + index: int, expected: int +) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(2) + condition = cirq.KeyCondition(cirq.MeasurementKey("control"), index=index) + result = simulator.run( + cirq.Circuit( + cirq.measure(qubits[0], key="control"), + cirq.X(qubits[0]), + cirq.measure(qubits[0], key="control"), + cirq.X(qubits[1]).with_classical_controls(condition), + cirq.measure(qubits[1], key="target"), + ) + ) + + np.testing.assert_array_equal(result.measurements["target"], [[expected]]) + assert simulator.read_register(qubits) == 1 + 2 * expected + + +@pytest.mark.parametrize("index", [1, -2]) +def test_key_condition_rejects_nonexistent_measurement_occurrence(index: int) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(2) + condition = cirq.KeyCondition(cirq.MeasurementKey("control"), index=index) + + with pytest.raises(IndexError): + simulator.run( + cirq.Circuit( + cirq.measure(qubits[0], key="control"), + cirq.X(qubits[1]).with_classical_controls(condition), + ) + ) + + +def test_classical_control_history_is_cleared_between_runs() -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(2) + simulator.run( + cirq.Circuit(cirq.X(qubits[0]), cirq.measure(qubits[0], key="control")) + ) + condition = cirq.KeyCondition(cirq.MeasurementKey("control"), index=0) + circuit = cirq.Circuit( + cirq.measure(qubits[0], key="control"), + cirq.X(qubits[1]).with_classical_controls(condition), + cirq.measure(qubits[1], key="target"), + ) + + result = _assert_matches_dense(simulator, circuit, qubits) + + np.testing.assert_array_equal(result.measurements["target"], [[0]]) + + +@pytest.mark.parametrize("first, second", list(itertools.product([0, 1], repeat=2))) +def test_multiple_classical_controls_are_conjoined(first: int, second: int) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(3) + circuit = cirq.Circuit( + _prepare_basis(qubits, first + 2 * second), + cirq.measure(qubits[0], key="a"), + cirq.measure(qubits[1], key="b"), + cirq.X(qubits[2]).with_classical_controls("a", "b"), + cirq.measure(qubits[2], key="target"), + ) + + result = _assert_matches_dense(simulator, circuit, qubits) + + np.testing.assert_array_equal(result.measurements["target"], [[first & second]]) + + +@pytest.mark.parametrize("initial", range(4)) +@pytest.mark.parametrize( + "expression", + [ + sympy.Eq(sympy.Symbol("control"), 1), + sympy.Eq(sympy.Symbol("control"), 2), + sympy.Lt(sympy.Symbol("control"), 2), + sympy.And( + sympy.Ne(sympy.Symbol("control"), 0), + sympy.Ne(sympy.Symbol("control"), 3), + ), + ], +) +def test_sympy_classical_conditions_use_cirq_integer_order( + initial: int, expression: sympy.Basic +) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(3) + circuit = cirq.Circuit( + _prepare_basis(qubits[:2], initial), + cirq.measure(*qubits[:2], key="control"), + cirq.X(qubits[2]).with_classical_controls(cirq.SympyCondition(expression)), + cirq.measure(qubits[2], key="target"), + ) + + _assert_matches_dense(simulator, circuit, qubits) + + +@pytest.mark.parametrize("initial", range(4)) +def test_sympy_condition_substitutes_multiple_measurement_keys(initial: int) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(3) + condition = cirq.SympyCondition( + sympy.Eq(sympy.Symbol("a") + sympy.Symbol("b"), 1) + ) + circuit = cirq.Circuit( + _prepare_basis(qubits[:2], initial), + cirq.measure(qubits[0], key="a"), + cirq.measure(qubits[1], key="b"), + cirq.X(qubits[2]).with_classical_controls(condition), + cirq.measure(qubits[2], key="target"), + ) + + _assert_matches_dense(simulator, circuit, qubits) + + +@pytest.mark.parametrize("invert_mask", [(True, False), (False, True)]) +def test_sympy_condition_uses_inverted_measurement( + invert_mask: tuple[bool, ...], +) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(3) + condition = cirq.SympyCondition(sympy.Eq(sympy.Symbol("control"), 2)) + circuit = cirq.Circuit( + cirq.measure(*qubits[:2], key="control", invert_mask=invert_mask), + cirq.X(qubits[2]).with_classical_controls(condition), + cirq.measure(qubits[2], key="target"), + ) + + _assert_matches_dense(simulator, circuit, qubits) + + + +@pytest.mark.parametrize("a, b", list(itertools.product(range(4), range(8)))) +def test_arithmetic_registers_are_little_endian_and_preserve_spectators( + a: int, b: int +) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(6) + first = [qubits[2], qubits[0]] + second = [qubits[4], qubits[1], qubits[3]] + gate = SparseSimArithmeticGate([2, 3], lambda x, y: (x, (x + y) % 8)) + assert cirq.num_qubits(gate) == 5 + result = simulator.run( + cirq.Circuit( + _prepare_basis(first, a), + _prepare_basis(second, b), + cirq.X(qubits[5]), + gate(*(first + second)), + cirq.measure(*first, key="a"), + cirq.measure(*second, key="b"), + cirq.measure(qubits[5], key="spectator"), + ) + ) + + expected_b = (a + b) % 8 + assert simulator.read_register(first) == a + assert simulator.read_register(second) == expected_b + assert simulator.read_register(qubits[5:]) == 1 + np.testing.assert_array_equal( + result.measurements["a"], [[(a >> i) & 1 for i in range(2)]] + ) + np.testing.assert_array_equal( + result.measurements["b"], [[(expected_b >> i) & 1 for i in range(3)]] + ) + np.testing.assert_array_equal(result.measurements["spectator"], [[1]]) + expected_bits = [ + (a >> 1) & 1, + (expected_b >> 1) & 1, + a & 1, + (expected_b >> 2) & 1, + expected_b & 1, + 1, + ] + np.testing.assert_allclose( + _state_vector(simulator, qubits), + np.eye(64)[cirq.big_endian_bits_to_int(expected_bits)], + ) + + +def test_arithmetic_preserves_complex_superposition() -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(4) + preparation = cirq.Circuit( + cirq.H.on_each(*qubits), cirq.T(qubits[0]), cirq.S(qubits[2]) + ) + gate = SparseSimArithmeticGate([2, 2], lambda a, b: (a, (a + b) % 4)) + permutation = np.zeros((16, 16), dtype=np.complex128) + for a, b in itertools.product(range(4), repeat=2): + permutation[4 * a + (a + b) % 4, 4 * a + b] = 1 + reference_circuit = preparation + cirq.Circuit( + cirq.MatrixGate(permutation)(qubits[1], qubits[0], qubits[3], qubits[2]) + ) + + simulator.run(preparation + cirq.Circuit(gate(*qubits))) + reference = cirq.Simulator(dtype=np.complex128).simulate( + reference_circuit, qubit_order=qubits + ) + + np.testing.assert_allclose( + _state_vector(simulator, qubits), reference.final_state_vector, atol=1e-12 + ) + assert len(simulator.basis_states) == 16 + + +def test_zero_width_arithmetic_registers() -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(2) + gate = SparseSimArithmeticGate( + [0, 2], lambda empty, value: (empty, (value + 1) % 4) + ) + zero_qubit_gate = SparseSimArithmeticGate([], lambda: ()) + + simulator.run( + cirq.Circuit(cirq.X.on_each(*qubits), gate(*qubits), zero_qubit_gate()) + ) + + assert cirq.num_qubits(zero_qubit_gate) == 0 + assert simulator.read_register(qubits) == 0 + np.testing.assert_allclose(_state_vector(simulator, qubits), [1, 0, 0, 0]) + + +def test_large_sparse_register_does_not_require_a_dense_state() -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(128) + gate = SparseSimArithmeticGate([128], lambda value: (value + 1,)) + expected = (1 << 127) + (1 << 63) + 1 + result = simulator.run( + cirq.Circuit( + cirq.X(qubits[127]), + cirq.CNOT(qubits[127], qubits[63]), + gate(*qubits), + cirq.measure(qubits[127], qubits[63], qubits[0], key="set_bits"), + ) + ) + + assert simulator.basis_states == [expected] + assert simulator.amplitudes == [1] + assert simulator.read_register(qubits) == expected + np.testing.assert_array_equal(result.measurements["set_bits"], [[1, 1, 1]]) + + +class _AncillaPhaseGate(cirq.Gate): + def __init__(self, *, release: bool = True, uncompute: bool = True) -> None: + self.release = release + self.uncompute = uncompute + + def _num_qubits_(self) -> int: + return 1 + + def _decompose_with_context_( + self, qubits: Sequence[cirq.Qid], context: cirq.DecompositionContext + ) -> Iterator[cirq.OP_TREE]: + ancilla = context.qubit_manager.qalloc(1)[0] + yield [cirq.CNOT(qubits[0], ancilla), [cirq.Z(ancilla)]] + if self.uncompute: + yield cirq.CNOT(qubits[0], ancilla) + if self.release: + context.qubit_manager.qfree([ancilla]) + + +class _NestedPhaseGate(cirq.Gate): + def _num_qubits_(self) -> int: + return 1 + + def _decompose_(self, qubits: Sequence[cirq.Qid]) -> Iterator[cirq.OP_TREE]: + yield [[_AncillaPhaseGate()(qubits[0])]] + + +def test_recursive_decomposition_allocates_uncomputes_and_reuses_ancilla() -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(1) + circuit = cirq.Circuit( + cirq.H(qubits[0]), + _NestedPhaseGate()(qubits[0]), + _NestedPhaseGate()(qubits[0]), + cirq.H(qubits[0]), + cirq.measure(qubits[0], key="m"), + ) + + result = simulator.run(circuit) + + np.testing.assert_array_equal(result.measurements["m"], [[0]]) + np.testing.assert_allclose(_state_vector(simulator, qubits), [1, 0], atol=1e-12) + assert simulator.qubit_manager.num_allocated_qubits() == 1 + assert simulator.qubit_manager.num_qubits == 2 + assert simulator.qubit_manager.qalloc(1) == [cirq.LineQubit(1)] + + +def test_tagged_and_circuit_operations() -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(2) + subcircuit = cirq.FrozenCircuit( + cirq.H(qubits[0]), cirq.CNOT(*qubits), cirq.S(qubits[1]) + ) + circuit = cirq.Circuit( + cirq.CircuitOperation(subcircuit).repeat(3), + cirq.X(qubits[0]).with_tags("tag"), + ) + + _assert_matches_dense(simulator, circuit, qubits) + + +def test_qubit_allocation_release_and_reuse() -> None: + simulator = SparseSimulator() + manager = simulator.qubit_manager + assert isinstance(manager, cirq.QubitManager) + assert manager.qalloc(0) == [] + assert manager.num_allocated_qubits() == 0 + qubits = manager.qalloc(4) + assert qubits == cirq.LineQubit.range(4) + assert manager.num_allocated_qubits() == 4 + + manager.qfree([qubits[1], qubits[3]]) + + assert manager.num_allocated_qubits() == 2 + assert manager.qalloc(3) == [qubits[3], qubits[1], cirq.LineQubit(4)] + assert manager.num_allocated_qubits() == 5 + assert manager.free_qubits == [] + assert manager.free_qubits_set == set() + manager.qfree([]) + assert manager.num_allocated_qubits() == 5 + + +def test_reused_qubits_simulate_with_noncontiguous_register_order() -> None: + simulator = SparseSimulator() + original = simulator.qubit_manager.qalloc(6) + simulator.qubit_manager.qfree([original[1], original[4]]) + recycled = simulator.qubit_manager.qalloc(2) + qubits = [recycled[1], original[5], recycled[0]] + circuit = cirq.Circuit( + cirq.H(qubits[0]), cirq.CNOT(qubits[0], qubits[2]), cirq.Y(qubits[1]) + ) + + _assert_matches_dense(simulator, circuit, qubits) + + +@pytest.mark.parametrize("gate", [cirq.X, cirq.H], ids=repr) +def test_qubits_cannot_be_freed_unless_every_branch_is_zero(gate: cirq.Gate) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(1) + simulator.run(cirq.Circuit(gate(*qubits))) + + with pytest.raises(RuntimeError, match="not in zero state"): + simulator.qubit_manager.qfree(qubits) + + assert simulator.qubit_manager.num_allocated_qubits() == 1 + assert simulator.qubit_manager.free_qubits_set == set() + + +def test_double_release_is_rejected() -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(1) + simulator.qubit_manager.qfree(qubits) + + with pytest.raises(RuntimeError, match="released twice"): + simulator.qubit_manager.qfree(qubits) + + assert simulator.qubit_manager.num_allocated_qubits() == 0 + assert simulator.qubit_manager.qalloc(2) == cirq.LineQubit.range(2) + + +@pytest.mark.parametrize( + "qubit", [cirq.NamedQubit("invalid"), cirq.LineQubit(-1), cirq.LineQubit(1)] +) +def test_free_rejects_invalid_qubits(qubit: cirq.Qid) -> None: + simulator = SparseSimulator() + simulator.qubit_manager.qalloc(1) + + with pytest.raises(AssertionError): + simulator.qubit_manager.qfree([qubit]) + + assert simulator.qubit_manager.num_allocated_qubits() == 1 + + +def test_nonbinary_allocation_and_borrowing_are_unsupported() -> None: + simulator = SparseSimulator() + + with pytest.raises(AssertionError): + simulator.qubit_manager.qalloc(1, dim=3) + with pytest.raises(NotImplementedError, match="qborrow is not supported"): + simulator.qubit_manager.qborrow(1) + + assert simulator.qubit_manager.num_allocated_qubits() == 0 + + +@pytest.mark.parametrize( + "gate, message", + [ + (_AncillaPhaseGate(release=False), "did not free all allocated qubits"), + (_AncillaPhaseGate(uncompute=False), "not in zero state"), + ], +) +def test_invalid_ancilla_decompositions_raise(gate: cirq.Gate, message: str) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(1) + + with pytest.raises(RuntimeError, match=message): + simulator.run(cirq.Circuit(cirq.X(qubits[0]), gate(*qubits))) + + +class _UnsupportedGate(cirq.Gate): + def _num_qubits_(self) -> int: + return 2 + + +@pytest.mark.parametrize( + "gate", + [ + _UnsupportedGate(), + cirq.depolarize(0.1), + cirq.amplitude_damp(0.1), + cirq.X ** sympy.Symbol("theta"), + cirq.CZ ** sympy.Symbol("theta"), + ], + ids=repr, +) +def test_unsupported_operations_raise_value_error(gate: cirq.Gate) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(cirq.num_qubits(gate)) + + with pytest.raises(ValueError, match="Operation cannot be simulated"): + simulator.run(cirq.Circuit(gate(*qubits))) + + +def test_unsupported_classical_condition_raises_value_error() -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(1) + condition = cirq.BitMaskKeyCondition(cirq.MeasurementKey("m"), bitmask=1) + circuit = cirq.Circuit( + cirq.measure(*qubits, key="m"), + cirq.X(*qubits).with_classical_controls(condition), + ) + + with pytest.raises(ValueError, match="Unsupported classical condition"): + simulator.run(circuit) + + +@pytest.mark.parametrize( + "condition", + ["missing", cirq.SympyCondition(sympy.Eq(sympy.Symbol("missing"), 1))], +) +def test_missing_classical_measurement_is_not_silently_ignored( + condition: str | cirq.Condition, +) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(1) + + with pytest.raises((KeyError, ValueError), match="missing"): + simulator.run( + cirq.Circuit(cirq.X(*qubits).with_classical_controls(condition)) + ) + + +@pytest.mark.parametrize("output", [(), (0, 0)]) +def test_arithmetic_gate_rejects_wrong_output_register_count( + output: tuple[int, ...], +) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(1) + gate = SparseSimArithmeticGate([1], lambda _: output) + + with pytest.raises(AssertionError): + simulator.run(cirq.Circuit(gate(*qubits))) \ No newline at end of file From 307f3a281790b6750ce57a7dbb73b8df209374bd Mon Sep 17 00:00:00 2001 From: Dmytro Fedoriaka Date: Sun, 27 Sep 2026 21:48:14 -0700 Subject: [PATCH 02/11] implement cirq.SparseSimulator --- cirq_sparse_sim/README.md | 45 ++- cirq_sparse_sim/sparse_sim.py | 450 +++++++++++++++++++++---- cirq_sparse_sim/sparse_sim_test.py | 512 ++++++++++++++++++++++++++++- 3 files changed, 927 insertions(+), 80 deletions(-) diff --git a/cirq_sparse_sim/README.md b/cirq_sparse_sim/README.md index bdb982b..0ab2d17 100644 --- a/cirq_sparse_sim/README.md +++ b/cirq_sparse_sim/README.md @@ -13,6 +13,39 @@ enough), so I just wrote this one for my needs. AI was used while writing this simulator. +## Usage + +```python +import cirq + +from cirq_sparse_sim.sparse_sim import SparseSimulator + +simulator = SparseSimulator(seed=42) +q0, q1 = simulator.qubit_manager.qalloc(2) +circuit = cirq.Circuit(cirq.H(q0), cirq.CNOT(q0, q1), cirq.measure(q0, q1, key="m")) + +samples = simulator.run(circuit, repetitions=100).measurements["m"] # (100, 2) +trial = simulator.simulate(circuit) +wavefunction = trial.final_state_vector +steps = list(simulator.simulate_moment_steps(circuit)) +``` + +The simulator implements `cirq.SimulatorBase` with sparse simulation states and +per-moment result snapshots. Each repetition starts in the all-zero state and +executes the same sparse gate and measurement logic. The simulator's +`basis_states`, `amplitudes`, `measurement_results`, and `read_register` expose +the last shot. Accessing a result's state vector explicitly materializes a dense +array; simulation itself remains sparse. Step-result sampling does not collapse +or otherwise change the saved state. + +Qubits must still be allocated through `qubit_manager`. Nonempty parameter +resolvers and non-integer initial states are not supported. Integer initial +states and returned wavefunctions use Cirq's big-endian `qubit_order`; +`read_register` and sparse basis indices retain their little-endian convention. +For backwards compatibility, `run` accepts circuits without measurements and +returns only the last occurrence of each measurement key in each shot, rather +than all repeated-key records. + ## Tests From the repository root, with Cirq, NumPy, SymPy, and pytest installed: @@ -22,12 +55,8 @@ python -m pytest cirq_sparse_sim/sparse_sim_test.py ``` The tests cover measurement results and final states, seeded random circuits -against Cirq's dense and Clifford simulators, arithmetic gates, classical -controls, reset, decomposition, qubit management, and error cases. Dense state +against Cirq's dense and Clifford simulators, independent repetitions, moment +steps and state snapshots, arithmetic gates, classical controls, reset, +decomposition, qubit management, and error cases. Dense state comparisons retain global phase; Clifford comparisons allow an overall global -phase because stabilizer representations need not preserve it. - -This suite targets the existing single-shot `run(circuit)` API and its -`cirq.Result` output. Full `cirq.SimulatorBase` compatibility, including -`simulate`, repetitions, parameter sweeps, and arbitrary qubit types, is not -implemented or claimed by these tests. \ No newline at end of file +phase because stabilizer representations need not preserve it. \ No newline at end of file diff --git a/cirq_sparse_sim/sparse_sim.py b/cirq_sparse_sim/sparse_sim.py index 484b8cb..3ec8491 100644 --- a/cirq_sparse_sim/sparse_sim.py +++ b/cirq_sparse_sim/sparse_sim.py @@ -1,4 +1,4 @@ -"""Sparse, single-shot simulator for a restricted subset of Cirq circuits. +"""Sparse SimulatorBase implementation for a restricted subset of Cirq circuits. The simulator keeps the quantum state as a sparse list of computational basis states and their amplitudes instead of a dense state vector. This is efficient @@ -10,15 +10,21 @@ 1. Start in ``|0...0>`` represented as one basis state with amplitude ``1``. 2. Apply operations one by one, updating sparse basis states/amplitudes. 3. Measurements collapse the sparse state stochastically according to Born rule. -4. Return a single-shot ``cirq.ResultDict``. +4. Repeat from the initial state for each requested shot. This simulator is intentionally limited and aimed at tests and debugging. Unsupported operations raise ``ValueError``. """ +from __future__ import annotations + +import copy import math +import operator import random -from collections.abc import Callable, Iterable, Sequence +from collections.abc import Callable, Iterable, Iterator, Mapping, Sequence +from types import NotImplementedType +from typing import Any, Self import cirq import numpy as np @@ -55,7 +61,7 @@ def num_qubits(self) -> int: class SparseSimQubitManager(cirq.QubitManager): """Simple qubit manager with allocation/release over ``LineQubit`` indices.""" - def __init__(self, simulator: "SparseSimulator") -> None: + def __init__(self, simulator: SparseSimulator | _SparseState) -> None: """Initialize the simulator.""" self.simulator = simulator self.num_qubits = 0 @@ -96,6 +102,11 @@ def num_allocated_qubits(self) -> int: """Returns number of allocated qubits.""" return self.num_qubits - len(self.free_qubits) + def _copy_allocation_from(self, other: SparseSimQubitManager) -> None: + self.num_qubits = other.num_qubits + self.free_qubits = other.free_qubits.copy() + self.free_qubits_set = other.free_qubits_set.copy() + def _apply_cx(state: int, control_bit: int, target_bit: int) -> int: if (state >> control_bit) % 2 == 1: @@ -120,44 +131,15 @@ def _normalize(amplitudes: list[complex]) -> list[complex]: return [a / norm for a in amplitudes] -class SparseSimulator: - """Sparse simulator for a restricted set of Cirq operations. - - Supported operations include: - - ``X``, ``CNOT``, ``CCNOT`` (basis permutation). - - ``CZPowGate`` with numeric exponent (phase on ``|11>`` branch). - - ``SparseSimArithmeticGate``. - - ``MeasurementGate`` and ``ResetChannel``. - - ``ClassicallyControlledOperation`` with ``KeyCondition`` or ``SympyCondition``. - - Any single-qubit gate with known unitary. - - Other operations are attempted via one-step decomposition and processed - recursively. If an operation still cannot be handled, simulation fails. - - Usage: - - Create a simulator instance: ``sim = SparseSimulator()``. - - Allocate qubits using ``sim.qubit_manager.qalloc``. Qubits created in any other - way are not supported. - - Create a circuit. - - Call ``sim.run(circuit)``. - - Use ``sim.read_register`` to read the value stored in a given register, - interpreted as an unsigned little-endian integer. - - Inspect measurement results in the ``cirq.Result`` object returned by - ``sim.run``. - - Notes: - - ``run`` is single-shot only. - - ``measurement_results`` stores the latest integer value per measurement key. - These values are little-endian; classical conditions use Cirq's big-endian - integers and measurement occurrence indices. - - ``read_register`` requires the final state to be a single basis state. - If it's not the case, add measurements. - - """ +class _SparseState(cirq.QuantumStateRepresentation, cirq.ClassicalDataStoreReader): + """The sparse execution engine and its classical measurement history.""" - def __init__(self, seed: int | None = None) -> None: - """Initialize SparseSimulator.""" - self.random = random.Random(seed) + def __init__( + self, + random_source: random.Random, + classical_data: cirq.ClassicalDataStoreReader | None = None, + ) -> None: + self.random: random.Random | np.random.RandomState = random_source self.qubit_manager = SparseSimQubitManager(self) # Sparse state is represented by a list of basis states with @@ -171,6 +153,23 @@ def __init__(self, seed: int | None = None) -> None: # Length of measurements (to restore measurements as bit vectors). self._meas_len: dict[str, int] = {} + self._records: dict[cirq.MeasurementKey, list[tuple[int, ...]]] = {} + if classical_data is not None: + if classical_data.channel_records: + raise NotImplementedError( + "Channel measurement records are not supported" + ) + for key, records in classical_data.records.items(): + self._records[key] = list(records) + if records: + name = str(key) + self._measurement_history[name] = [ + cirq.big_endian_bits_to_int(bits) for bits in records + ] + self._meas_len[name] = len(records[-1]) + self.measurement_results[name] = sum( + bit << i for i, bit in enumerate(records[-1]) + ) def _apply_x(self, qid: int) -> None: self.basis_states = [s ^ (1 << qid) for s in self.basis_states] @@ -313,6 +312,9 @@ def _apply_measurement_gate( self._measurement_history.setdefault(key, []).append( cirq.big_endian_bits_to_int(measured_bits) ) + self._records.setdefault(cirq.MeasurementKey.parse_serialized(key), []).append( + tuple(measured_bits) + ) def _apply_reset(self, target_qubit: int) -> None: """Reset one qubit to ``|0>`` via measure-and-conditional-flip.""" @@ -398,39 +400,347 @@ def _run_op(self, op: cirq.Operation, context: cirq.DecompositionContext) -> Non raise ValueError(f"Operation cannot be simulated: {op}.") - def run(self, circuit: cirq.Circuit) -> cirq.Result: - """Run ``circuit`` once and return ``cirq.ResultDict``. - - Returned ``measurements[key]`` is a boolean array with shape - ``(1, num_measured_qubits_for_key)``. - """ - self.basis_states = [0] - self.amplitudes = [1.0 + 0.0j] - self.measurement_results = {} - self._measurement_history = {} - self._meas_len = {} - # Circuit must use only qubits allocated using qubit manager of this simulator. - context = cirq.DecompositionContext(self.qubit_manager) - for op in circuit.all_operations(): - self._run_op(op, context) - - measurements = { - key: np.asarray( - [[((value >> i) & 1) == 1 for i in range(self._meas_len[key])]], - dtype=np.bool_, - ) - for key, value in self.measurement_results.items() - } - return cirq.ResultDict(params=cirq.ParamResolver({}), measurements=measurements) - def read_register(self, register: Sequence[cirq.Qid]) -> int: """Read a little-endian register as an integer.""" - assert len(self.basis_states) == 1, ( - "Final state is superposition, add measurements." - ) + assert ( + len(self.basis_states) == 1 + ), "Final state is superposition, add measurements." s = self.basis_states[0] return sum(((s >> q.x) & 1) << i for i, q in enumerate(register)) def _is_qubit_zero(self, qid: int) -> bool: """Checks whether qubit is in 0 state.""" - return all((s >> qid) & 1 == 0 for s in self.basis_states) \ No newline at end of file + return all((s >> qid) & 1 == 0 for s in self.basis_states) + + def copy(self, deep_copy_buffers: bool = True) -> Self: + """Copy the state; there are no reusable scratch buffers to share.""" + result = copy.copy(self) + result.basis_states = self.basis_states.copy() + result.amplitudes = self.amplitudes.copy() + result.measurement_results = self.measurement_results.copy() + result._measurement_history = { + key: values.copy() for key, values in self._measurement_history.items() + } + result._meas_len = self._meas_len.copy() + result._records = {key: values.copy() for key, values in self._records.items()} + result.qubit_manager = SparseSimQubitManager(result) + result.qubit_manager._copy_allocation_from(self.qubit_manager) + return result + + def measure( + self, axes: Sequence[int], seed: cirq.RANDOM_STATE_OR_SEED_LIKE = None + ) -> list[int]: + if seed is not None: + self.random = cirq.value.parse_random_state(seed) + return [self._measure_qubit(axis) for axis in axes] + + def sample( + self, + axes: Sequence[int], + repetitions: int = 1, + seed: cirq.RANDOM_STATE_OR_SEED_LIKE = None, + ) -> np.ndarray: + repetitions = operator.index(repetitions) + if repetitions < 0: + raise ValueError("repetitions must be non-negative") + return super().sample(axes, repetitions, seed).reshape(repetitions, len(axes)) + + @property + def records(self) -> Mapping[cirq.MeasurementKey, list[tuple[int, ...]]]: + return self._records + + @property + def channel_records(self) -> Mapping[cirq.MeasurementKey, list[int]]: + return {} + + def keys(self) -> tuple[cirq.MeasurementKey, ...]: + return tuple(self._records) + + def get_digits(self, key: cirq.MeasurementKey, index: int = -1) -> tuple[int, ...]: + return self._records[key][index] + + def get_int(self, key: cirq.MeasurementKey, index: int = -1) -> int: + return cirq.big_endian_bits_to_int(self.get_digits(key, index)) + + +class _SparseSimulationState(cirq.SimulationState[_SparseState]): + """Cirq's simulation-state interface over the existing sparse engine.""" + + @property + def sparse_state(self) -> _SparseState: + return self._state + + @property + def classical_data(self) -> cirq.ClassicalDataStoreReader: + return self._state + + def get_axes(self, qubits: Sequence[cirq.Qid]) -> list[int]: + axes = [] + for qubit in qubits: + if qubit not in self.qubit_map: + raise ValueError(f"Qubit {qubit} is not in this simulation state") + if not isinstance(qubit, cirq.LineQubit): + raise ValueError("Only LineQubits are supported") + axes.append(qubit.x) + return axes + + def apply_operation(self, op: cirq.Operation) -> None: + self.get_axes(op.qubits) + self._state._run_op(op, cirq.DecompositionContext(self._state.qubit_manager)) + + def _act_on_fallback_( + self, action: Any, qubits: Sequence[cirq.Qid], allow_decompose: bool = True + ) -> bool | NotImplementedType: + if isinstance(action, cirq.Gate): + action = action.on(*qubits) + if not isinstance(action, cirq.Operation): + return NotImplemented + self.apply_operation(action) + return True + + def measure( + self, + qubits: Sequence[cirq.Qid], + key: str, + invert_mask: Sequence[bool], + confusion_map: dict[tuple[int, ...], np.ndarray], + ) -> None: + self.apply_operation( + cirq.measure( + *qubits, + key=key, + invert_mask=tuple(invert_mask), + confusion_map=confusion_map, + ) + ) + + def state_vector(self) -> np.ndarray: + axes = self.get_axes(self.qubits) + vector = np.zeros(1 << len(axes), dtype=np.complex128) + included_bits = sum(1 << axis for axis in axes) + for basis, amplitude in zip( + self._state.basis_states, self._state.amplitudes, strict=True + ): + if basis & ~included_bits: + raise ValueError("An omitted qubit is not in the zero state") + index = 0 + for axis in axes: + index = (index << 1) | ((basis >> axis) & 1) + vector[index] = amplitude + return vector + + +class SparseSimulatorStep(cirq.StepResultBase[_SparseSimulationState]): + """A snapshot after a moment, with non-collapsing sampling support.""" + + def state_vector(self, copy: bool = True) -> np.ndarray: + """Materialize the sparse state in the simulation's qubit order.""" + return self._merged_sim_state.state_vector() + + +class SparseSimulatorTrialResult( + cirq.SimulationTrialResultBase[_SparseSimulationState] +): + """Final sparse simulation state, materialized as a vector only on request.""" + + @property + def final_state_vector(self) -> np.ndarray: + return self._get_merged_sim_state().state_vector() + + +class SparseSimulator( + cirq.SimulatorBase[ + SparseSimulatorStep, SparseSimulatorTrialResult, _SparseSimulationState + ] +): + """Sparse simulator for a restricted set of Cirq operations. + + Supports X, CNOT, CCNOT, numeric CZ powers, single-qubit unitaries, + SparseSimArithmeticGate, measurement, reset, and classical controls using + KeyCondition or SympyCondition. Other operations are tried by decomposition. + + Allocate qubits through ``qubit_manager.qalloc`` before constructing circuits. + ``run(circuit, repetitions=n)`` executes independent shots starting in zero. + ``simulate`` and ``simulate_moment_steps`` expose final and intermediate + states, with integer initial states interpreted in Cirq's qubit order. + Nonempty parameter resolvers are not supported. + + ``basis_states``, ``amplitudes``, and ``measurement_results`` describe the + latest shot. Register reads and measurement_results are little-endian; + Cirq state vectors and classical conditions use big-endian ordering. + ``read_register`` requires a single basis state, so measure first if needed. + + For backwards compatibility, run also accepts circuits without measurements + and returns only the latest occurrence of each measurement key per shot. + """ + + def __init__(self, seed: int | None = None) -> None: + super().__init__(dtype=np.complex128, seed=seed, split_untangled_states=False) + self.random = random.Random(seed) + self._state = _SparseState(self.random) + self.qubit_manager = SparseSimQubitManager(self) + + @property + def basis_states(self) -> list[int]: + return self._state.basis_states + + @property + def amplitudes(self) -> list[complex]: + return self._state.amplitudes + + @property + def measurement_results(self) -> dict[str, int]: + return self._state.measurement_results + + def read_register(self, register: Sequence[cirq.Qid]) -> int: + """Read a little-endian register from the latest shot.""" + return self._state.read_register(register) + + def _is_qubit_zero(self, qid: int) -> bool: + return self._state._is_qubit_zero(qid) + + def _create_partial_simulation_state( + self, + initial_state: Any, + qubits: Sequence[cirq.Qid], + classical_data: cirq.ClassicalDataStore, + ) -> _SparseSimulationState: + if not isinstance(initial_state, (int, np.integer)): + raise NotImplementedError("Only integer initial states are supported") + if not 0 <= initial_state < (1 << len(qubits)): + raise ValueError("initial_state is out of range for the supplied qubits") + basis = 0 + for i, qubit in enumerate(qubits): + if ( + not isinstance(qubit, cirq.LineQubit) + or not 0 <= qubit.x < self.qubit_manager.num_qubits + or qubit.x in self.qubit_manager.free_qubits_set + ): + raise ValueError("Circuit qubits must be allocated by qubit_manager") + basis |= ((int(initial_state) >> (len(qubits) - i - 1)) & 1) << qubit.x + state = _SparseState(self.random, classical_data) + state.basis_states = [basis] + state.qubit_manager._copy_allocation_from(self.qubit_manager) + self._state = state + return _SparseSimulationState( + state=state, qubits=qubits, classical_data=classical_data, prng=self._prng + ) + + def _create_step_result( + self, sim_state: cirq.SimulationStateBase[_SparseSimulationState] + ) -> SparseSimulatorStep: + self._state = sim_state.create_merged_state().sparse_state + self.qubit_manager._copy_allocation_from(self._state.qubit_manager) + return SparseSimulatorStep(sim_state.copy()) + + def _create_simulator_trial_result( + self, + params: cirq.ParamResolver, + measurements: dict[str, np.ndarray], + final_simulator_state: cirq.SimulationStateBase[_SparseSimulationState], + ) -> SparseSimulatorTrialResult: + return SparseSimulatorTrialResult(params, measurements, final_simulator_state) + + def _can_be_in_run_prefix(self, val: Any) -> bool: + # Preserve per-shot execution and the original decomposition context. + return False + + def _core_iterator( + self, + circuit: cirq.AbstractCircuit, + sim_state: cirq.SimulationStateBase[_SparseSimulationState], + all_measurements_are_terminal: bool = False, + ) -> Iterator[SparseSimulatorStep]: + state = sim_state.create_merged_state() + # Keep the engine's dispatch, including its supported classical controls. + for moment in circuit if len(circuit) else [cirq.Moment()]: + self._state = state.sparse_state + try: + for op in moment.operations: + if all_measurements_are_terminal and cirq.is_measurement(op): + continue + state.apply_operation(op) + finally: + self.qubit_manager._copy_allocation_from(self._state.qubit_manager) + yield self._create_step_result(state) + + @staticmethod + def _check_parameters( + param_resolver: cirq.ParamResolverOrSimilarType, + ) -> None: + if cirq.ParamResolver(param_resolver): + raise NotImplementedError("Parameter resolvers are not supported") + + def simulate_moment_steps( + self, + circuit: cirq.AbstractCircuit, + param_resolver: cirq.ParamResolverOrSimilarType = None, + qubit_order: cirq.QubitOrderOrList = cirq.QubitOrder.DEFAULT, + initial_state: Any = None, + ) -> Iterator[SparseSimulatorStep]: + self._check_parameters(param_resolver) + return super().simulate_moment_steps( + circuit, param_resolver, qubit_order, initial_state + ) + + def simulate_sweep_iter( + self, + program: cirq.AbstractCircuit, + params: cirq.Sweepable, + qubit_order: cirq.QubitOrderOrList = cirq.QubitOrder.DEFAULT, + initial_state: Any = None, + ) -> Iterator[SparseSimulatorTrialResult]: + resolvers = list(cirq.to_resolvers(params)) + for resolver in resolvers: + self._check_parameters(resolver) + return super().simulate_sweep_iter( + program, resolvers, qubit_order, initial_state + ) + + def run_sweep_iter( + self, + program: cirq.AbstractCircuit, + params: cirq.Sweepable, + repetitions: int = 1, + ) -> Iterator[cirq.Result]: + for resolver in cirq.to_resolvers(params): + yield cirq.ResultDict( + params=resolver, + records=self._run(program, resolver, repetitions), + ) + + def _run( + self, + circuit: cirq.AbstractCircuit, + param_resolver: cirq.ParamResolver, + repetitions: int, + ) -> dict[str, np.ndarray]: + self._check_parameters(param_resolver) + repetitions = operator.index(repetitions) + if repetitions < 0: + raise ValueError("repetitions must be non-negative") + if repetitions == 0: + operations = cirq.decompose( + circuit, + keep=lambda op: isinstance(op.gate, cirq.MeasurementGate) + or not cirq.is_measurement(op), + ) + return { + cirq.measurement_key_name(op): np.empty( + (0, 1, len(op.qubits)), dtype=np.bool_ + ) + for op in operations + if isinstance(op.gate, cirq.MeasurementGate) + } + + records: dict[str, list[Sequence[int]]] = {} + qubits = tuple(sorted(circuit.all_qubits())) + for _ in range(repetitions): + for step in self._base_iterator(circuit, qubits, 0): + pass + for key, bits in step.measurements.items(): + records.setdefault(key, []).append(bits) + return { + key: np.asarray(values, dtype=np.bool_)[:, np.newaxis, :] + for key, values in records.items() + } diff --git a/cirq_sparse_sim/sparse_sim_test.py b/cirq_sparse_sim/sparse_sim_test.py index 3ad5f30..ef2343b 100644 --- a/cirq_sparse_sim/sparse_sim_test.py +++ b/cirq_sparse_sim/sparse_sim_test.py @@ -1,4 +1,4 @@ -"""Tests for SparseSimulator's existing single-shot API.""" +"""Tests for sparse simulation, Cirq interfaces, and independent repetitions.""" import itertools import math @@ -1068,4 +1068,512 @@ def test_arithmetic_gate_rejects_wrong_output_register_count( gate = SparseSimArithmeticGate([1], lambda _: output) with pytest.raises(AssertionError): - simulator.run(cirq.Circuit(gate(*qubits))) \ No newline at end of file + simulator.run(cirq.Circuit(gate(*qubits))) + + +@pytest.mark.parametrize("repetitions", [0, 1, 2, 17]) +def test_simulator_base_run_supports_repetitions(repetitions: int) -> None: + simulator = SparseSimulator(seed=1) + qubits = simulator.qubit_manager.qalloc(3) + circuit = cirq.Circuit( + cirq.X(qubits[0]), + cirq.measure(qubits[2], qubits[0], qubits[1], key="bits"), + cirq.measure(qubits[0], key="flag"), + ) + assert isinstance(simulator, cirq.SimulatorBase) + sampler: cirq.Sampler = simulator + + result = sampler.run(program=circuit, repetitions=repetitions) + + assert result.params == cirq.ParamResolver({}) + assert result.measurements["bits"].dtype == np.bool_ + assert result.records["bits"].shape == (repetitions, 1, 3) + np.testing.assert_array_equal( + result.measurements["bits"], np.tile([0, 1, 0], (repetitions, 1)) + ) + np.testing.assert_array_equal( + result.measurements["flag"], np.ones((repetitions, 1), dtype=bool) + ) + if repetitions: + assert simulator.read_register(qubits) == 1 + assert simulator.measurement_results == {"bits": 2, "flag": 1} + + +@pytest.mark.parametrize( + "kind", ["terminal", "mid-circuit", "reset", "arithmetic", "ancilla"] +) +def test_repetitions_match_successive_single_shots(kind: str) -> None: + batched = SparseSimulator(seed=3456) + sequential = SparseSimulator(seed=3456) + qubits = batched.qubit_manager.qalloc(3) + assert sequential.qubit_manager.qalloc(3) == qubits + circuit = cirq.Circuit( + cirq.ry(0.71)(qubits[0]), cirq.H(qubits[1]), cirq.CNOT(*qubits[1:]) + ) + if kind == "mid-circuit": + circuit += cirq.Circuit( + cirq.measure(qubits[0], key="control"), + cirq.X(qubits[2]).with_classical_controls("control"), + cirq.H(qubits[0]), + ) + elif kind == "reset": + circuit += cirq.Circuit(cirq.reset(qubits[1]), cirq.H(qubits[2])) + elif kind == "arithmetic": + gate = SparseSimArithmeticGate([3], lambda value: ((value + 3) % 8,)) + circuit += cirq.Circuit(gate(*qubits)) + elif kind == "ancilla": + circuit += cirq.Circuit(_NestedPhaseGate()(qubits[0])) + circuit += cirq.measure(*qubits, key="m") + repetitions = 64 + manager = batched.qubit_manager + + result = batched.run(circuit, repetitions=repetitions) + individual = [sequential.run(circuit) for _ in range(repetitions)] + + assert result.measurements.keys() == individual[0].measurements.keys() + for key in result.measurements: + np.testing.assert_array_equal( + result.measurements[key], + np.concatenate([shot.measurements[key] for shot in individual]), + ) + assert batched.basis_states == sequential.basis_states + np.testing.assert_allclose(batched.amplitudes, sequential.amplitudes) + assert batched.measurement_results == sequential.measurement_results + assert len({tuple(row) for row in result.measurements["m"]}) > 1 + assert batched.qubit_manager is manager + assert manager.num_allocated_qubits() == 3 + assert manager.num_qubits == (4 if kind == "ancilla" else 3) + + +def test_repetitions_preserve_entangled_measurement_correlations() -> None: + simulator = SparseSimulator(seed=5) + qubits = simulator.qubit_manager.qalloc(2) + circuit = cirq.Circuit( + cirq.H(qubits[0]), + cirq.CNOT(*qubits), + cirq.measure(qubits[0], key="a"), + cirq.measure(qubits[1], key="b"), + ) + + result = simulator.run(circuit, repetitions=100) + + np.testing.assert_array_equal(result.measurements["a"], result.measurements["b"]) + assert set(result.measurements["a"].ravel()) == {False, True} + assert simulator.read_register(qubits) == 3 * int(result.measurements["a"][-1, 0]) + + +def test_repetitions_clear_classical_history_between_shots() -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(2) + first_measurement = cirq.KeyCondition(cirq.MeasurementKey("control"), index=0) + circuit = cirq.Circuit( + cirq.measure(qubits[0], key="control"), + cirq.X(qubits[0]), + cirq.measure(qubits[0], key="control"), + cirq.X(qubits[1]).with_classical_controls(first_measurement), + cirq.measure(qubits[1], key="target"), + ) + + result = simulator.run(circuit, repetitions=7) + + np.testing.assert_array_equal(result.measurements["control"], np.ones((7, 1))) + np.testing.assert_array_equal(result.measurements["target"], np.zeros((7, 1))) + assert simulator.read_register(qubits) == 1 + + +def test_repetitions_retain_latest_repeated_key_even_when_width_changes() -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(2) + circuit = cirq.Circuit( + cirq.measure(qubits[0], key="m"), + cirq.X(qubits[1]), + cirq.measure(*qubits, key="m"), + ) + + result = simulator.run(circuit, repetitions=5) + + assert result.records["m"].shape == (5, 1, 2) + np.testing.assert_array_equal(result.measurements["m"], np.tile([0, 1], (5, 1))) + + +def test_zero_repetitions_preserve_state_results_and_random_stream() -> None: + simulator = SparseSimulator(seed=5) + qubits = simulator.qubit_manager.qalloc(2) + simulator.run(cirq.Circuit(cirq.X(qubits[0]), cirq.measure(*qubits, key="old"))) + random_state = simulator.random.getstate() + circuit = cirq.Circuit( + cirq.H(qubits[0]), + cirq.CircuitOperation(cirq.FrozenCircuit(cirq.measure(*qubits, key="new"))), + ) + + result = simulator.run(circuit, repetitions=0) + + assert result.measurements["new"].shape == (0, 2) + assert simulator.basis_states == [1] + assert simulator.measurement_results == {"old": 1} + assert simulator.random.getstate() == random_state + + +@pytest.mark.parametrize("repetitions", [1, 3]) +def test_repetitions_without_measurements_preserve_legacy_run( + repetitions: int, +) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(1) + + result = simulator.run(cirq.Circuit(cirq.X(*qubits)), repetitions=repetitions) + + assert result.measurements == {} + assert simulator.read_register(qubits) == 1 + + +@pytest.mark.parametrize("repetitions", [-1, -3]) +def test_negative_repetitions_are_rejected(repetitions: int) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(1) + + with pytest.raises(ValueError, match="repetitions"): + simulator.run(cirq.Circuit(cirq.measure(*qubits)), repetitions=repetitions) + + +def test_parameterless_sweeps_use_repetitions() -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(1) + circuit = cirq.Circuit(cirq.X(*qubits), cirq.measure(*qubits, key="m")) + + results = simulator.run_sweep(circuit, params=[{}, {}], repetitions=4) + + assert len(results) == 2 + for result in results: + np.testing.assert_array_equal(result.measurements["m"], np.ones((4, 1))) + + +@pytest.mark.parametrize("initial_state", [0, 1, 3, 5, np.int64(6)]) +@pytest.mark.parametrize("reverse_order", [False, True]) +def test_simulate_respects_qubit_order_and_integer_initial_state( + initial_state: int, reverse_order: bool +) -> None: + simulator = SparseSimulator(seed=0) + allocated = simulator.qubit_manager.qalloc(6) + qubits = [allocated[4], allocated[0], allocated[2]] + order = qubits[::-1] if reverse_order else qubits + circuit = cirq.Circuit( + cirq.ry(0.4)(qubits[0]), cirq.CNOT(qubits[0], qubits[1]), cirq.T(qubits[2]) + ) + + actual = simulator.simulate(circuit, qubit_order=order, initial_state=initial_state) + expected = cirq.Simulator(dtype=np.complex128).simulate( + circuit, qubit_order=order, initial_state=int(initial_state) + ) + + assert isinstance(actual, cirq.SimulationTrialResultBase) + assert actual.params == cirq.ParamResolver({}) + assert actual.measurements == {} + assert actual.qubit_map == {q: i for i, q in enumerate(order)} + np.testing.assert_allclose(actual.final_state_vector, expected.final_state_vector) + np.testing.assert_allclose( + actual.final_state_vector, _state_vector(simulator, order) + ) + + +@pytest.mark.parametrize("seed", range(5)) +def test_simulate_random_circuits_matches_dense_simulator(seed: int) -> None: + simulator = SparseSimulator(seed=seed) + qubits = simulator.qubit_manager.qalloc(4) + circuit = cirq.testing.random_circuit( + qubits, n_moments=20, op_density=0.8, random_state=seed + ) + + result = simulator.simulate(circuit, qubit_order=qubits) + expected = cirq.Simulator(dtype=np.complex128).simulate(circuit, qubit_order=qubits) + + np.testing.assert_allclose( + result.final_state_vector, expected.final_state_vector, atol=1e-8 + ) + + +def test_simulate_measurements_and_final_state_are_independent_snapshots() -> None: + simulator = SparseSimulator(seed=1) + qubits = simulator.qubit_manager.qalloc(2) + result = simulator.simulate( + cirq.Circuit( + cirq.H(qubits[0]), + cirq.CNOT(*qubits), + cirq.measure(*qubits, key="m"), + ) + ) + measured = result.measurements["m"].copy() + state_vector = result.final_state_vector.copy() + + simulator.run(cirq.Circuit(cirq.X(qubits[0])), repetitions=2) + + assert measured.shape == (2,) + assert measured[0] == measured[1] + np.testing.assert_array_equal(result.measurements["m"], measured) + np.testing.assert_allclose(result.final_state_vector, state_vector) + np.testing.assert_allclose(state_vector, np.eye(4)[3 * measured[0]]) + + +def test_moment_steps_match_dense_states_and_do_not_alias() -> None: + simulator = SparseSimulator(seed=0) + qubits = simulator.qubit_manager.qalloc(2) + circuit = cirq.Circuit( + cirq.Moment(cirq.H(qubits[0])), + cirq.Moment(cirq.T(qubits[0])), + cirq.Moment(cirq.CNOT(*qubits)), + cirq.Moment(cirq.measure(qubits[0], key="m")), + ) + + steps = list(simulator.simulate_moment_steps(circuit, qubit_order=qubits)) + + assert len(steps) == len(circuit) + assert all(isinstance(step, cirq.StepResultBase) for step in steps) + for i in range(3): + expected = cirq.Simulator(dtype=np.complex128).simulate( + circuit[: i + 1], qubit_order=qubits + ) + assert steps[i].measurements == {} + np.testing.assert_allclose(steps[i].state_vector(), expected.final_state_vector) + bit = steps[-1].measurements["m"][0] + expected_final = np.eye(4, dtype=np.complex128)[3 * bit] + expected_final *= np.exp(1j * np.pi * bit / 4) + np.testing.assert_allclose(steps[-1].state_vector(), expected_final) + first_vector = steps[0].state_vector() + first_vector[:] = 0 + assert np.linalg.norm(steps[0].state_vector()) == pytest.approx(1) + + +@pytest.mark.parametrize("repetitions", [0, 1, 50]) +def test_step_sampling_is_seeded_and_does_not_collapse_state(repetitions: int) -> None: + simulator = SparseSimulator(seed=0) + qubits = simulator.qubit_manager.qalloc(3) + circuit = cirq.Circuit( + cirq.X(qubits[2]), cirq.H(qubits[0]), cirq.CNOT(*qubits[:2]) + ) + step = list(simulator.simulate_moment_steps(circuit))[-1] + before = step.state_vector() + random_state = simulator.random.getstate() + + samples = step.sample([qubits[2], qubits[1], qubits[0]], repetitions, seed=42) + + assert samples.shape == (repetitions, 3) + np.testing.assert_array_equal(samples[:, 0], np.ones(repetitions)) + np.testing.assert_array_equal(samples[:, 1], samples[:, 2]) + np.testing.assert_array_equal( + samples, step.sample([qubits[2], qubits[1], qubits[0]], repetitions, seed=42) + ) + np.testing.assert_allclose(step.state_vector(), before) + assert step.measurements == {} + assert simulator.random.getstate() == random_state + assert step.sample([], repetitions, seed=42).shape == (repetitions, 0) + with pytest.raises(ValueError, match="repetitions"): + step.sample(qubits, repetitions=-1) + + +def test_simulation_state_copy_is_independent_and_supports_act_on() -> None: + simulator = SparseSimulator(seed=0) + qubits = simulator.qubit_manager.qalloc(2) + result = simulator.simulate( + cirq.Circuit(cirq.H(qubits[0]), cirq.measure(qubits[1], key="old")) + ) + state = result.get_state_containing_qubit(qubits[0]) + copied = state.copy(deep_copy_buffers=False) + before = state.state_vector() + + cirq.act_on(cirq.X(qubits[1]), copied) + cirq.act_on(cirq.measure(qubits[1], key="new"), copied) + + assert copied.classical_data.get_int(cirq.MeasurementKey("new")) == 1 + assert copied.classical_data.records[cirq.MeasurementKey("new")] == [(1,)] + assert state.classical_data.keys() == (cirq.MeasurementKey("old"),) + assert copied.classical_data.channel_records == {} + assert state.sparse_state.qubit_manager is not copied.sparse_state.qubit_manager + np.testing.assert_allclose(state.state_vector(), before) + np.testing.assert_allclose(result.final_state_vector, before) + cirq.act_on(cirq.X, copied, qubits=[qubits[1]]) + np.testing.assert_allclose(copied.state_vector(), before) + with pytest.raises(ValueError, match="not in this simulation state"): + cirq.act_on(cirq.X(cirq.LineQubit(100)), copied) + with pytest.raises(TypeError, match="Failed to act"): + cirq.act_on(object(), copied, qubits=qubits) + + +def test_empty_simulation_and_moment_steps_return_initial_state() -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(2) + + result = simulator.simulate(cirq.Circuit(), qubit_order=qubits, initial_state=2) + steps = list( + simulator.simulate_moment_steps( + cirq.Circuit(), qubit_order=qubits, initial_state=2 + ) + ) + + assert len(steps) == 1 + assert result.measurements == {} + assert steps[0].measurements == {} + np.testing.assert_array_equal(result.final_state_vector, [0, 0, 1, 0]) + np.testing.assert_array_equal(steps[0].state_vector(), result.final_state_vector) + + +def test_simulate_keeps_large_register_sparse_until_vector_is_requested() -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(128) + + result = simulator.simulate( + cirq.Circuit(cirq.X(qubits[127])), qubit_order=qubits + ) + + state = result.get_state_containing_qubit(qubits[127]).sparse_state + assert state.basis_states == [1 << 127] + assert state.amplitudes == [1] + assert simulator.read_register(qubits) == 1 << 127 + + +@pytest.mark.parametrize("method", ["run", "simulate", "simulate_moment_steps"]) +@pytest.mark.parametrize("resolver", [None, {}, cirq.ParamResolver({})]) +def test_empty_parameter_resolvers_are_accepted(method: str, resolver) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(1) + circuit = cirq.Circuit(cirq.X(*qubits), cirq.measure(*qubits, key="m")) + + result = getattr(simulator, method)(circuit, param_resolver=resolver) + if method == "simulate_moment_steps": + result = list(result)[-1] + + assert np.all(result.measurements["m"]) + + +@pytest.mark.parametrize( + "method", + ["run", "simulate", "simulate_moment_steps", "run_sweep", "simulate_sweep"], +) +@pytest.mark.parametrize( + "resolver", [{"theta": 0.5}, cirq.ParamResolver({"theta": 0.5})] +) +def test_nonempty_parameter_resolvers_are_rejected(method: str, resolver) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(1) + circuit = cirq.Circuit( + (cirq.X ** sympy.Symbol("theta"))(*qubits), + cirq.measure(*qubits, key="m"), + ) + parameter = "params" if "sweep" in method else "param_resolver" + + with pytest.raises(NotImplementedError, match="Parameter resolvers"): + result = getattr(simulator, method)(circuit, **{parameter: resolver}) + if method == "simulate_moment_steps": + list(result) + + +@pytest.mark.parametrize("initial_state", [-1, 4]) +def test_out_of_range_initial_states_are_rejected(initial_state: int) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(2) + + with pytest.raises(ValueError, match="initial_state"): + simulator.simulate( + cirq.Circuit(), qubit_order=qubits, initial_state=initial_state + ) + + +def test_dense_initial_states_are_explicitly_unsupported() -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(1) + + with pytest.raises(NotImplementedError, match="integer initial states"): + simulator.simulate( + cirq.Circuit(), qubit_order=qubits, initial_state=np.array([1, 0]) + ) + + +def test_failed_decomposition_keeps_public_allocator_in_sync() -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(1) + + with pytest.raises(RuntimeError, match="did not free all allocated qubits"): + simulator.run( + cirq.Circuit(cirq.X(*qubits), _AncillaPhaseGate(release=False)(*qubits)) + ) + + assert simulator.qubit_manager.num_allocated_qubits() == 2 + assert simulator.qubit_manager.qalloc(1) == [cirq.LineQubit(2)] + + +def test_parameterless_simulation_sweep_has_independent_results() -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(1) + circuit = cirq.Circuit(cirq.X(*qubits), cirq.measure(*qubits, key="m")) + + results = simulator.simulate_sweep(circuit, params=[{}, {}]) + + assert len(results) == 2 + for result in results: + np.testing.assert_array_equal(result.measurements["m"], [1]) + np.testing.assert_allclose(result.final_state_vector, [0, 1]) + assert ( + results[0].get_state_containing_qubit(qubits[0]) + is not results[1].get_state_containing_qubit(qubits[0]) + ) + + +def test_simulation_can_continue_from_returned_state() -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(1) + first = simulator.simulate(cirq.Circuit(cirq.H(*qubits))) + state = first.get_state_containing_qubit(qubits[0]) + + second = simulator.simulate( + cirq.Circuit(cirq.H(*qubits)), initial_state=state.copy() + ) + + np.testing.assert_allclose(second.final_state_vector, [1, 0], atol=1e-12) + np.testing.assert_allclose(first.final_state_vector, [2**-0.5, 2**-0.5]) + + +def test_unallocated_circuit_qubits_are_rejected() -> None: + simulator = SparseSimulator() + simulator.qubit_manager.qalloc(1) + + with pytest.raises(ValueError, match="allocated by qubit_manager"): + simulator.run(cirq.Circuit(cirq.X(cirq.LineQubit(1)))) + + +def test_partial_simulation_state_preserves_supplied_classical_data() -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(2) + data = cirq.ClassicalDataDictionaryStore() + key = cirq.MeasurementKey("control") + data.record_measurement(key, [0, 1], qubits) + state = simulator._create_partial_simulation_state(0, qubits, data) + + cirq.act_on(cirq.X(qubits[0]).with_classical_controls(str(key)), state) + + assert isinstance(state, cirq.SimulationState) + assert state.classical_data.get_int(key) == 1 + assert state.classical_data.get_digits(key) == (0, 1) + assert state.sparse_state.measurement_results["control"] == 2 + np.testing.assert_allclose(state.state_vector(), [0, 0, 1, 0]) + + +def test_partial_state_rejects_unsupported_channel_records() -> None: + simulator = SparseSimulator() + data = cirq.ClassicalDataDictionaryStore() + data.record_channel_measurement(cirq.MeasurementKey("channel"), 1) + + with pytest.raises(NotImplementedError, match="Channel measurement"): + simulator._create_partial_simulation_state(0, (), data) + + +@pytest.mark.parametrize("method", ["run", "simulate"]) +def test_rejected_parameters_do_not_clear_previous_state(method: str) -> None: + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(1) + circuit = cirq.Circuit(cirq.X(*qubits), cirq.measure(*qubits, key="m")) + simulator.run(circuit) + + with pytest.raises(NotImplementedError, match="Parameter resolvers"): + getattr(simulator, method)(circuit, param_resolver={"unused": 1}) + + assert simulator.basis_states == [1] + assert simulator.measurement_results == {"m": 1} \ No newline at end of file From d8612d5b4d21365e40d1aa7ec17ed45c876cf453 Mon Sep 17 00:00:00 2001 From: Dmytro Fedoriaka Date: Sun, 27 Sep 2026 21:52:11 -0700 Subject: [PATCH 03/11] support params --- cirq_sparse_sim/sparse_sim.py | 26 ++++------------------- cirq_sparse_sim/sparse_sim_test.py | 34 ++++++++++++++---------------- 2 files changed, 20 insertions(+), 40 deletions(-) diff --git a/cirq_sparse_sim/sparse_sim.py b/cirq_sparse_sim/sparse_sim.py index 3ec8491..46500a9 100644 --- a/cirq_sparse_sim/sparse_sim.py +++ b/cirq_sparse_sim/sparse_sim.py @@ -29,6 +29,7 @@ import cirq import numpy as np import sympy +from cirq.sim.simulator import check_all_resolved _AMPLITUDE_EPS = 1e-9 @@ -563,7 +564,6 @@ class SparseSimulator( ``run(circuit, repetitions=n)`` executes independent shots starting in zero. ``simulate`` and ``simulate_moment_steps`` expose final and intermediate states, with integer initial states interpreted in Cirq's qubit order. - Nonempty parameter resolvers are not supported. ``basis_states``, ``amplitudes``, and ``measurement_results`` describe the latest shot. Register reads and measurement_results are little-endian; @@ -664,25 +664,6 @@ def _core_iterator( self.qubit_manager._copy_allocation_from(self._state.qubit_manager) yield self._create_step_result(state) - @staticmethod - def _check_parameters( - param_resolver: cirq.ParamResolverOrSimilarType, - ) -> None: - if cirq.ParamResolver(param_resolver): - raise NotImplementedError("Parameter resolvers are not supported") - - def simulate_moment_steps( - self, - circuit: cirq.AbstractCircuit, - param_resolver: cirq.ParamResolverOrSimilarType = None, - qubit_order: cirq.QubitOrderOrList = cirq.QubitOrder.DEFAULT, - initial_state: Any = None, - ) -> Iterator[SparseSimulatorStep]: - self._check_parameters(param_resolver) - return super().simulate_moment_steps( - circuit, param_resolver, qubit_order, initial_state - ) - def simulate_sweep_iter( self, program: cirq.AbstractCircuit, @@ -692,7 +673,7 @@ def simulate_sweep_iter( ) -> Iterator[SparseSimulatorTrialResult]: resolvers = list(cirq.to_resolvers(params)) for resolver in resolvers: - self._check_parameters(resolver) + check_all_resolved(cirq.resolve_parameters(program, resolver)) return super().simulate_sweep_iter( program, resolvers, qubit_order, initial_state ) @@ -715,7 +696,8 @@ def _run( param_resolver: cirq.ParamResolver, repetitions: int, ) -> dict[str, np.ndarray]: - self._check_parameters(param_resolver) + circuit = cirq.resolve_parameters(circuit, param_resolver) + check_all_resolved(circuit) repetitions = operator.index(repetitions) if repetitions < 0: raise ValueError("repetitions must be non-negative") diff --git a/cirq_sparse_sim/sparse_sim_test.py b/cirq_sparse_sim/sparse_sim_test.py index ef2343b..cfc8713 100644 --- a/cirq_sparse_sim/sparse_sim_test.py +++ b/cirq_sparse_sim/sparse_sim_test.py @@ -1017,8 +1017,6 @@ def _num_qubits_(self) -> int: _UnsupportedGate(), cirq.depolarize(0.1), cirq.amplitude_damp(0.1), - cirq.X ** sympy.Symbol("theta"), - cirq.CZ ** sympy.Symbol("theta"), ], ids=repr, ) @@ -1444,26 +1442,22 @@ def test_empty_parameter_resolvers_are_accepted(method: str, resolver) -> None: assert np.all(result.measurements["m"]) -@pytest.mark.parametrize( - "method", - ["run", "simulate", "simulate_moment_steps", "run_sweep", "simulate_sweep"], -) -@pytest.mark.parametrize( - "resolver", [{"theta": 0.5}, cirq.ParamResolver({"theta": 0.5})] -) -def test_nonempty_parameter_resolvers_are_rejected(method: str, resolver) -> None: +def test_parameter_sweeps_resolve_gates() -> None: simulator = SparseSimulator() qubits = simulator.qubit_manager.qalloc(1) circuit = cirq.Circuit( (cirq.X ** sympy.Symbol("theta"))(*qubits), cirq.measure(*qubits, key="m"), ) - parameter = "params" if "sweep" in method else "param_resolver" + params = [{"theta": 0}, {"theta": 1}] - with pytest.raises(NotImplementedError, match="Parameter resolvers"): - result = getattr(simulator, method)(circuit, **{parameter: resolver}) - if method == "simulate_moment_steps": - list(result) + run_results = simulator.run_sweep(circuit, params=params, repetitions=2) + simulation_results = simulator.simulate_sweep(circuit, params=params) + + np.testing.assert_array_equal(run_results[0].measurements["m"], [[0], [0]]) + np.testing.assert_array_equal(run_results[1].measurements["m"], [[1], [1]]) + np.testing.assert_array_equal(simulation_results[0].measurements["m"], [0]) + np.testing.assert_array_equal(simulation_results[1].measurements["m"], [1]) @pytest.mark.parametrize("initial_state", [-1, 4]) @@ -1566,14 +1560,18 @@ def test_partial_state_rejects_unsupported_channel_records() -> None: @pytest.mark.parametrize("method", ["run", "simulate"]) -def test_rejected_parameters_do_not_clear_previous_state(method: str) -> None: +def test_unresolved_parameters_do_not_clear_previous_state(method: str) -> None: simulator = SparseSimulator() qubits = simulator.qubit_manager.qalloc(1) circuit = cirq.Circuit(cirq.X(*qubits), cirq.measure(*qubits, key="m")) simulator.run(circuit) + parameterized_circuit = cirq.Circuit( + (cirq.X ** sympy.Symbol("theta"))(*qubits), + cirq.measure(*qubits, key="m"), + ) - with pytest.raises(NotImplementedError, match="Parameter resolvers"): - getattr(simulator, method)(circuit, param_resolver={"unused": 1}) + with pytest.raises(ValueError, match="symbols were not specified"): + getattr(simulator, method)(parameterized_circuit) assert simulator.basis_states == [1] assert simulator.measurement_results == {"m": 1} \ No newline at end of file From 320f6376a52246d96632757b5724254fddfdce55 Mon Sep 17 00:00:00 2001 From: Dmytro Fedoriaka Date: Sun, 27 Sep 2026 21:57:11 -0700 Subject: [PATCH 04/11] support multiple measurments per key --- cirq_sparse_sim/sparse_sim.py | 42 ++++++++++++++++-------------- cirq_sparse_sim/sparse_sim_test.py | 29 ++++++++++++--------- 2 files changed, 40 insertions(+), 31 deletions(-) diff --git a/cirq_sparse_sim/sparse_sim.py b/cirq_sparse_sim/sparse_sim.py index 46500a9..524f77c 100644 --- a/cirq_sparse_sim/sparse_sim.py +++ b/cirq_sparse_sim/sparse_sim.py @@ -570,8 +570,7 @@ class SparseSimulator( Cirq state vectors and classical conditions use big-endian ordering. ``read_register`` requires a single basis state, so measure first if needed. - For backwards compatibility, run also accepts circuits without measurements - and returns only the latest occurrence of each measurement key per shot. + Accepts circuits without measurements. """ def __init__(self, seed: int | None = None) -> None: @@ -701,28 +700,33 @@ def _run( repetitions = operator.index(repetitions) if repetitions < 0: raise ValueError("repetitions must be non-negative") + operations = cirq.decompose( + circuit, + keep=lambda op: isinstance(op.gate, cirq.MeasurementGate) + or not cirq.is_measurement(op), + ) + measurement_widths: dict[str, list[int]] = {} + for op in operations: + if isinstance(op.gate, cirq.MeasurementGate): + key = cirq.measurement_key_name(op) + measurement_widths.setdefault(key, []).append(len(op.qubits)) + for key, widths in measurement_widths.items(): + if len(set(widths)) != 1: + raise ValueError( + f"Different qid shapes for repeated measurement: key={key!r}" + ) + if repetitions == 0: - operations = cirq.decompose( - circuit, - keep=lambda op: isinstance(op.gate, cirq.MeasurementGate) - or not cirq.is_measurement(op), - ) return { - cirq.measurement_key_name(op): np.empty( - (0, 1, len(op.qubits)), dtype=np.bool_ - ) - for op in operations - if isinstance(op.gate, cirq.MeasurementGate) + key: np.empty((0, len(widths), widths[0]), dtype=np.bool_) + for key, widths in measurement_widths.items() } - records: dict[str, list[Sequence[int]]] = {} + records: dict[str, list[Sequence[Sequence[int]]]] = {} qubits = tuple(sorted(circuit.all_qubits())) for _ in range(repetitions): for step in self._base_iterator(circuit, qubits, 0): pass - for key, bits in step.measurements.items(): - records.setdefault(key, []).append(bits) - return { - key: np.asarray(values, dtype=np.bool_)[:, np.newaxis, :] - for key, values in records.items() - } + for key, values in step._classical_data.records.items(): + records.setdefault(str(key), []).append(values) + return {key: np.asarray(values, dtype=np.bool_) for key, values in records.items()} diff --git a/cirq_sparse_sim/sparse_sim_test.py b/cirq_sparse_sim/sparse_sim_test.py index cfc8713..74d087d 100644 --- a/cirq_sparse_sim/sparse_sim_test.py +++ b/cirq_sparse_sim/sparse_sim_test.py @@ -451,21 +451,26 @@ def test_measurement_keys_and_result_interface( np.testing.assert_array_equal(_state_vector(simulator, qubits), np.eye(8)[4]) -def test_repeated_measurement_key_retains_latest_value() -> None: +def test_repeated_measurement_key_retains_all_occurrences() -> None: simulator = SparseSimulator() qubits = simulator.qubit_manager.qalloc(2) result = simulator.run( cirq.Circuit( cirq.X(qubits[0]), - cirq.measure(qubits[0], key="m"), + cirq.measure(*qubits, key="m"), cirq.X(qubits[0]), cirq.X(qubits[1]), cirq.measure(*qubits, key="m"), - ) + ), + repetitions=2, ) assert simulator.measurement_results == {"m": 2} - np.testing.assert_array_equal(result.measurements["m"], [[0, 1]]) + np.testing.assert_array_equal( + result.records["m"], [[[1, 0], [0, 1]], [[1, 0], [0, 1]]] + ) + with pytest.raises(ValueError, match="repeated keys"): + _ = result.measurements assert simulator.read_register(qubits) == 2 @@ -621,7 +626,7 @@ def test_key_condition_index_selects_measurement_occurrence( ) ) - np.testing.assert_array_equal(result.measurements["target"], [[expected]]) + np.testing.assert_array_equal(result.records["target"], [[[expected]]]) assert simulator.read_register(qubits) == 1 + 2 * expected @@ -1174,12 +1179,14 @@ def test_repetitions_clear_classical_history_between_shots() -> None: result = simulator.run(circuit, repetitions=7) - np.testing.assert_array_equal(result.measurements["control"], np.ones((7, 1))) - np.testing.assert_array_equal(result.measurements["target"], np.zeros((7, 1))) + np.testing.assert_array_equal( + result.records["control"], np.tile([[[0], [1]]], (7, 1, 1)) + ) + np.testing.assert_array_equal(result.records["target"], np.zeros((7, 1, 1))) assert simulator.read_register(qubits) == 1 -def test_repetitions_retain_latest_repeated_key_even_when_width_changes() -> None: +def test_repeated_measurement_key_rejects_different_widths() -> None: simulator = SparseSimulator() qubits = simulator.qubit_manager.qalloc(2) circuit = cirq.Circuit( @@ -1188,10 +1195,8 @@ def test_repetitions_retain_latest_repeated_key_even_when_width_changes() -> Non cirq.measure(*qubits, key="m"), ) - result = simulator.run(circuit, repetitions=5) - - assert result.records["m"].shape == (5, 1, 2) - np.testing.assert_array_equal(result.measurements["m"], np.tile([0, 1], (5, 1))) + with pytest.raises(ValueError, match="Different qid shapes"): + simulator.run(circuit, repetitions=5) def test_zero_repetitions_preserve_state_results_and_random_stream() -> None: From f6ab3f5143aa20584ee7ecec0b7c76a3956a07d1 Mon Sep 17 00:00:00 2001 From: Dmytro Fedoriaka Date: Sun, 27 Sep 2026 22:30:21 -0700 Subject: [PATCH 05/11] support external qubits --- cirq_sparse_sim/sparse_sim.py | 121 +++++++++++++++++++++-------- cirq_sparse_sim/sparse_sim_test.py | 68 ++++++++++++++-- 2 files changed, 150 insertions(+), 39 deletions(-) diff --git a/cirq_sparse_sim/sparse_sim.py b/cirq_sparse_sim/sparse_sim.py index 524f77c..1644587 100644 --- a/cirq_sparse_sim/sparse_sim.py +++ b/cirq_sparse_sim/sparse_sim.py @@ -60,7 +60,7 @@ def num_qubits(self) -> int: class SparseSimQubitManager(cirq.QubitManager): - """Simple qubit manager with allocation/release over ``LineQubit`` indices.""" + """Qubit manager used for legacy allocation and decomposition ancillas.""" def __init__(self, simulator: SparseSimulator | _SparseState) -> None: """Initialize the simulator.""" @@ -70,18 +70,30 @@ def __init__(self, simulator: SparseSimulator | _SparseState) -> None: self.free_qubits_set: set[int] = set() def _allocate_qubit(self) -> cirq.LineQubit: - if len(self.free_qubits) == 0: + for index in range(len(self.free_qubits) - 1, -1, -1): + qubit_id = self.free_qubits[index] + qubit = cirq.LineQubit(qubit_id) + if not self.simulator._has_qubit(qubit): + del self.free_qubits[index] + self.free_qubits_set.remove(qubit_id) + return qubit + + while True: qubit_id = self.num_qubits self.num_qubits += 1 - else: - qubit_id = self.free_qubits.pop() - self.free_qubits_set.remove(qubit_id) - return cirq.LineQubit(qubit_id) + qubit = cirq.LineQubit(qubit_id) + if not self.simulator._has_qubit(qubit): + return qubit + self.free_qubits.append(qubit_id) + self.free_qubits_set.add(qubit_id) def qalloc(self, n: int, dim: int = 2) -> list[cirq.Qid]: """Allocate qubits.""" assert dim == 2 - return [self._allocate_qubit() for _ in range(n)] + qubits = [self._allocate_qubit() for _ in range(n)] + for qubit in qubits: + self.simulator._add_qubit(qubit) + return qubits def qborrow(self, n: int, dim: int = 2) -> list[cirq.Qid]: """Not implemented.""" @@ -94,8 +106,9 @@ def qfree(self, qubits: Iterable[cirq.Qid]) -> None: assert 0 <= q.x < self.num_qubits if q.x in self.free_qubits_set: raise RuntimeError(f"Qubit {q.x} released twice") - if not self.simulator._is_qubit_zero(q.x): + if not self.simulator._is_qubit_zero(q): raise RuntimeError(f"Qubit {q.x} released not in zero state") + self.simulator._remove_qubit(q) self.free_qubits.append(q.x) self.free_qubits_set.add(q.x) @@ -103,6 +116,11 @@ def num_allocated_qubits(self) -> int: """Returns number of allocated qubits.""" return self.num_qubits - len(self.free_qubits) + def allocated_qubits(self) -> Iterator[cirq.LineQubit]: + for qubit_id in range(self.num_qubits): + if qubit_id not in self.free_qubits_set: + yield cirq.LineQubit(qubit_id) + def _copy_allocation_from(self, other: SparseSimQubitManager) -> None: self.num_qubits = other.num_qubits self.free_qubits = other.free_qubits.copy() @@ -138,9 +156,13 @@ class _SparseState(cirq.QuantumStateRepresentation, cirq.ClassicalDataStoreReade def __init__( self, random_source: random.Random, + qubits: Sequence[cirq.Qid] = (), classical_data: cirq.ClassicalDataStoreReader | None = None, ) -> None: self.random: random.Random | np.random.RandomState = random_source + self.axis_by_qubit = {qubit: axis for axis, qubit in enumerate(qubits)} + self._next_axis = len(qubits) + self._free_axes: list[int] = [] self.qubit_manager = SparseSimQubitManager(self) # Sparse state is represented by a list of basis states with @@ -172,6 +194,20 @@ def __init__( bit << i for i, bit in enumerate(records[-1]) ) + def _has_qubit(self, qubit: cirq.Qid) -> bool: + return qubit in self.axis_by_qubit + + def _add_qubit(self, qubit: cirq.Qid) -> None: + if qubit in self.axis_by_qubit: + raise RuntimeError(f"Qubit {qubit} is already in the simulation state") + axis = self._free_axes.pop() if self._free_axes else self._next_axis + if axis == self._next_axis: + self._next_axis += 1 + self.axis_by_qubit[qubit] = axis + + def _remove_qubit(self, qubit: cirq.Qid) -> None: + self._free_axes.append(self.axis_by_qubit.pop(qubit)) + def _apply_x(self, qid: int) -> None: self.basis_states = [s ^ (1 << qid) for s in self.basis_states] @@ -350,7 +386,7 @@ def _run_op(self, op: cirq.Operation, context: cirq.DecompositionContext) -> Non return op_gate = op.gate - qubit_ids = [q.x for q in op.qubits] + qubit_ids = [self.axis_by_qubit[q] for q in op.qubits] if isinstance(op_gate, cirq.XPowGate) and op_gate.exponent == 1: self._apply_x(qubit_ids[0]) @@ -407,11 +443,20 @@ def read_register(self, register: Sequence[cirq.Qid]) -> int: len(self.basis_states) == 1 ), "Final state is superposition, add measurements." s = self.basis_states[0] - return sum(((s >> q.x) & 1) << i for i, q in enumerate(register)) + try: + return sum( + ((s >> self.axis_by_qubit[q]) & 1) << i + for i, q in enumerate(register) + ) + except KeyError as ex: + raise ValueError( + f"Qubit {ex.args[0]} is not in this simulation state" + ) from ex - def _is_qubit_zero(self, qid: int) -> bool: + def _is_qubit_zero(self, qubit: cirq.Qid) -> bool: """Checks whether qubit is in 0 state.""" - return all((s >> qid) & 1 == 0 for s in self.basis_states) + axis = self.axis_by_qubit.get(qubit) + return axis is None or all((s >> axis) & 1 == 0 for s in self.basis_states) def copy(self, deep_copy_buffers: bool = True) -> Self: """Copy the state; there are no reusable scratch buffers to share.""" @@ -424,6 +469,9 @@ def copy(self, deep_copy_buffers: bool = True) -> Self: } result._meas_len = self._meas_len.copy() result._records = {key: values.copy() for key, values in self._records.items()} + result.axis_by_qubit = self.axis_by_qubit.copy() + result._next_axis = self._next_axis + result._free_axes = self._free_axes.copy() result.qubit_manager = SparseSimQubitManager(result) result.qubit_manager._copy_allocation_from(self.qubit_manager) return result @@ -476,14 +524,12 @@ def classical_data(self) -> cirq.ClassicalDataStoreReader: return self._state def get_axes(self, qubits: Sequence[cirq.Qid]) -> list[int]: - axes = [] - for qubit in qubits: - if qubit not in self.qubit_map: - raise ValueError(f"Qubit {qubit} is not in this simulation state") - if not isinstance(qubit, cirq.LineQubit): - raise ValueError("Only LineQubits are supported") - axes.append(qubit.x) - return axes + try: + return [self._state.axis_by_qubit[qubit] for qubit in qubits] + except KeyError as ex: + raise ValueError( + f"Qubit {ex.args[0]} is not in this simulation state" + ) from ex def apply_operation(self, op: cirq.Operation) -> None: self.get_axes(op.qubits) @@ -560,7 +606,6 @@ class SparseSimulator( SparseSimArithmeticGate, measurement, reset, and classical controls using KeyCondition or SympyCondition. Other operations are tried by decomposition. - Allocate qubits through ``qubit_manager.qalloc`` before constructing circuits. ``run(circuit, repetitions=n)`` executes independent shots starting in zero. ``simulate`` and ``simulate_moment_steps`` expose final and intermediate states, with integer initial states interpreted in Cirq's qubit order. @@ -595,8 +640,17 @@ def read_register(self, register: Sequence[cirq.Qid]) -> int: """Read a little-endian register from the latest shot.""" return self._state.read_register(register) - def _is_qubit_zero(self, qid: int) -> bool: - return self._state._is_qubit_zero(qid) + def _has_qubit(self, qubit: cirq.Qid) -> bool: + return self._state._has_qubit(qubit) + + def _add_qubit(self, qubit: cirq.Qid) -> None: + self._state._add_qubit(qubit) + + def _remove_qubit(self, qubit: cirq.Qid) -> None: + self._state._remove_qubit(qubit) + + def _is_qubit_zero(self, qubit: cirq.Qid) -> bool: + return self._state._is_qubit_zero(qubit) def _create_partial_simulation_state( self, @@ -608,16 +662,21 @@ def _create_partial_simulation_state( raise NotImplementedError("Only integer initial states are supported") if not 0 <= initial_state < (1 << len(qubits)): raise ValueError("initial_state is out of range for the supplied qubits") + for qubit in qubits: + if qubit.dimension != 2: + raise ValueError("Only dimension-2 qubits are supported") + state_qubits = list(qubits) + circuit_qubits = set(qubits) + state_qubits.extend( + qubit + for qubit in self.qubit_manager.allocated_qubits() + if qubit not in circuit_qubits + ) + state = _SparseState(self.random, state_qubits, classical_data) basis = 0 for i, qubit in enumerate(qubits): - if ( - not isinstance(qubit, cirq.LineQubit) - or not 0 <= qubit.x < self.qubit_manager.num_qubits - or qubit.x in self.qubit_manager.free_qubits_set - ): - raise ValueError("Circuit qubits must be allocated by qubit_manager") - basis |= ((int(initial_state) >> (len(qubits) - i - 1)) & 1) << qubit.x - state = _SparseState(self.random, classical_data) + axis = state.axis_by_qubit[qubit] + basis |= ((int(initial_state) >> (len(qubits) - i - 1)) & 1) << axis state.basis_states = [basis] state.qubit_manager._copy_allocation_from(self.qubit_manager) self._state = state diff --git a/cirq_sparse_sim/sparse_sim_test.py b/cirq_sparse_sim/sparse_sim_test.py index 74d087d..a9809a7 100644 --- a/cirq_sparse_sim/sparse_sim_test.py +++ b/cirq_sparse_sim/sparse_sim_test.py @@ -16,10 +16,7 @@ def _state_vector( simulator: SparseSimulator, qubit_order: Sequence[cirq.Qid] ) -> np.ndarray: """Convert little-endian sparse indices to Cirq's specified tensor order.""" - qubit_ids = [] - for qubit in qubit_order: - assert isinstance(qubit, cirq.LineQubit) - qubit_ids.append(qubit.x) + qubit_ids = [simulator._state.axis_by_qubit[qubit] for qubit in qubit_order] assert len(simulator.basis_states) == len(simulator.amplitudes) assert len(set(simulator.basis_states)) == len(simulator.basis_states) @@ -899,6 +896,34 @@ def test_recursive_decomposition_allocates_uncomputes_and_reuses_ancilla() -> No assert simulator.qubit_manager.qalloc(1) == [cirq.LineQubit(1)] +def test_recursive_decomposition_supports_external_qubits() -> None: + simulator = SparseSimulator() + qubit = cirq.NamedQubit("external") + circuit = cirq.Circuit( + cirq.H(qubit), + _NestedPhaseGate()(qubit), + _NestedPhaseGate()(qubit), + cirq.H(qubit), + cirq.measure(qubit, key="m"), + ) + + result = simulator.run(circuit) + + np.testing.assert_array_equal(result.measurements["m"], [[0]]) + assert simulator._state.axis_by_qubit == {qubit: 0} + assert simulator.qubit_manager.num_allocated_qubits() == 0 + + +def test_internal_ancilla_does_not_alias_external_line_qubit() -> None: + simulator = SparseSimulator() + qubit = cirq.LineQubit(0) + + simulator.run(cirq.Circuit(_NestedPhaseGate()(qubit))) + + assert simulator._state.axis_by_qubit == {qubit: 0} + assert simulator.qubit_manager.qalloc(1) == [cirq.LineQubit(1)] + + def test_tagged_and_circuit_operations() -> None: simulator = SparseSimulator() qubits = simulator.qubit_manager.qalloc(2) @@ -1530,12 +1555,39 @@ def test_simulation_can_continue_from_returned_state() -> None: np.testing.assert_allclose(first.final_state_vector, [2**-0.5, 2**-0.5]) -def test_unallocated_circuit_qubits_are_rejected() -> None: +@pytest.mark.parametrize( + "qubits", + [ + cirq.NamedQubit.range(2, prefix="q"), + [cirq.GridQubit(4, 7), cirq.GridQubit(-2, 10)], + [cirq.LineQubit(100), cirq.LineQubit(10_000)], + ], +) +def test_arbitrary_circuit_qubits_use_compact_axes( + qubits: Sequence[cirq.Qid], +) -> None: + simulator = SparseSimulator() + circuit = cirq.Circuit( + cirq.X(qubits[0]), + cirq.CNOT(*qubits), + cirq.measure(*qubits, key="m"), + ) + + result = simulator.run(circuit) + + np.testing.assert_array_equal(result.measurements["m"], [[1, 1]]) + assert simulator._state.axis_by_qubit == { + qubit: axis for axis, qubit in enumerate(sorted(qubits)) + } + assert simulator.basis_states == [3] + + +def test_nonbinary_circuit_qubits_are_rejected() -> None: simulator = SparseSimulator() - simulator.qubit_manager.qalloc(1) + qubit = cirq.NamedQid("qutrit", dimension=3) - with pytest.raises(ValueError, match="allocated by qubit_manager"): - simulator.run(cirq.Circuit(cirq.X(cirq.LineQubit(1)))) + with pytest.raises(ValueError, match="dimension-2"): + simulator.run(cirq.Circuit(cirq.XPowGate(dimension=3)(qubit))) def test_partial_simulation_state_preserves_supplied_classical_data() -> None: From cfde06b15cd0034aa617efbe16b511b0c25fcd41 Mon Sep 17 00:00:00 2001 From: Dmytro Fedoriaka Date: Sun, 27 Sep 2026 23:02:01 -0700 Subject: [PATCH 06/11] support external qubits --- cirq_sparse_sim/README.md | 17 ++++++++-------- cirq_sparse_sim/sparse_sim.py | 38 ++++++++++++++++++++++++++--------- 2 files changed, 37 insertions(+), 18 deletions(-) diff --git a/cirq_sparse_sim/README.md b/cirq_sparse_sim/README.md index 0ab2d17..beff4cf 100644 --- a/cirq_sparse_sim/README.md +++ b/cirq_sparse_sim/README.md @@ -21,7 +21,7 @@ import cirq from cirq_sparse_sim.sparse_sim import SparseSimulator simulator = SparseSimulator(seed=42) -q0, q1 = simulator.qubit_manager.qalloc(2) +q0, q1 = cirq.NamedQubit.range(2, prefix="q") circuit = cirq.Circuit(cirq.H(q0), cirq.CNOT(q0, q1), cirq.measure(q0, q1, key="m")) samples = simulator.run(circuit, repetitions=100).measurements["m"] # (100, 2) @@ -38,13 +38,14 @@ the last shot. Accessing a result's state vector explicitly materializes a dense array; simulation itself remains sparse. Step-result sampling does not collapse or otherwise change the saved state. -Qubits must still be allocated through `qubit_manager`. Nonempty parameter -resolvers and non-integer initial states are not supported. Integer initial -states and returned wavefunctions use Cirq's big-endian `qubit_order`; -`read_register` and sparse basis indices retain their little-endian convention. -For backwards compatibility, `run` accepts circuits without measurements and -returns only the last occurrence of each measurement key in each shot, rather -than all repeated-key records. +Circuits may use arbitrary dimension-2 Cirq Qids; `qubit_manager` remains +available for backwards compatibility and decomposition ancillas. Integer +initial states and returned wavefunctions use Cirq's big-endian `qubit_order`. +`read_register` is little-endian in the supplied register order, while sparse +basis-state bit positions use the simulation's compact internal qubit mapping. +Parameter resolvers and repeated measurement-key records follow Cirq's normal +contracts. Non-integer initial states are not supported. For backwards +compatibility, `run` accepts circuits without measurements. ## Tests diff --git a/cirq_sparse_sim/sparse_sim.py b/cirq_sparse_sim/sparse_sim.py index 1644587..3a2c811 100644 --- a/cirq_sparse_sim/sparse_sim.py +++ b/cirq_sparse_sim/sparse_sim.py @@ -69,10 +69,23 @@ def __init__(self, simulator: SparseSimulator | _SparseState) -> None: self.free_qubits: list[int] = [] self.free_qubits_set: set[int] = set() - def _allocate_qubit(self) -> cirq.LineQubit: + def _qid_for_id(self, qubit_id: int) -> cirq.Qid: + if isinstance(self.simulator, _SparseState): + return cirq.ops.CleanQubit(qubit_id, prefix="_sparse_sim") + return cirq.LineQubit(qubit_id) + + @staticmethod + def _qubit_id(qubit: cirq.Qid) -> int: + if isinstance(qubit, cirq.LineQubit): + return qubit.x + if isinstance(qubit, cirq.ops.CleanQubit): + return qubit.id + raise AssertionError(f"Qubit {qubit} was not allocated by this manager") + + def _allocate_qubit(self) -> cirq.Qid: for index in range(len(self.free_qubits) - 1, -1, -1): qubit_id = self.free_qubits[index] - qubit = cirq.LineQubit(qubit_id) + qubit = self._qid_for_id(qubit_id) if not self.simulator._has_qubit(qubit): del self.free_qubits[index] self.free_qubits_set.remove(qubit_id) @@ -81,7 +94,7 @@ def _allocate_qubit(self) -> cirq.LineQubit: while True: qubit_id = self.num_qubits self.num_qubits += 1 - qubit = cirq.LineQubit(qubit_id) + qubit = self._qid_for_id(qubit_id) if not self.simulator._has_qubit(qubit): return qubit self.free_qubits.append(qubit_id) @@ -102,15 +115,17 @@ def qborrow(self, n: int, dim: int = 2) -> list[cirq.Qid]: def qfree(self, qubits: Iterable[cirq.Qid]) -> None: """Free qubits.""" for q in qubits: - assert isinstance(q, cirq.LineQubit) - assert 0 <= q.x < self.num_qubits - if q.x in self.free_qubits_set: - raise RuntimeError(f"Qubit {q.x} released twice") + qubit_id = self._qubit_id(q) + assert 0 <= qubit_id < self.num_qubits + if qubit_id in self.free_qubits_set: + raise RuntimeError(f"Qubit {qubit_id} released twice") + if not self.simulator._has_qubit(q): + raise RuntimeError(f"Qubit {q} is not in the simulation state") if not self.simulator._is_qubit_zero(q): - raise RuntimeError(f"Qubit {q.x} released not in zero state") + raise RuntimeError(f"Qubit {qubit_id} released not in zero state") self.simulator._remove_qubit(q) - self.free_qubits.append(q.x) - self.free_qubits_set.add(q.x) + self.free_qubits.append(qubit_id) + self.free_qubits_set.add(qubit_id) def num_allocated_qubits(self) -> int: """Returns number of allocated qubits.""" @@ -524,6 +539,9 @@ def classical_data(self) -> cirq.ClassicalDataStoreReader: return self._state def get_axes(self, qubits: Sequence[cirq.Qid]) -> list[int]: + for qubit in qubits: + if qubit not in self.qubit_map: + raise ValueError(f"Qubit {qubit} is not in this simulation state") try: return [self._state.axis_by_qubit[qubit] for qubit in qubits] except KeyError as ex: From b27599ebb756dcc6915c46776761870a956f0207 Mon Sep 17 00:00:00 2001 From: Dmytro Fedoriaka Date: Sun, 27 Sep 2026 23:14:40 -0700 Subject: [PATCH 07/11] add demo --- README.md | 2 + cirq_sparse_sim/Demo.ipynb | 231 +++++++++++++++++++++++++++++++++++++ 2 files changed, 233 insertions(+) create mode 100644 cirq_sparse_sim/Demo.ipynb diff --git a/README.md b/README.md index 2367cf4..78ec3e9 100644 --- a/README.md +++ b/README.md @@ -19,6 +19,8 @@ List of projects: * [Quantum Lookups](quantum_lookups/Quantum%20Lookups.ipynb) (2025) * [Screening Task for QOSF Mentorship Program](qosf_tasks/2025-gate-tomography/Gate%20Tomography.ipynb) (2025) * [Formal verification of Q# circuits](qsharp-verification/formal-verification-qsharp.ipynb) (2026) +* [Cirq Sparse Simulator](cirq_sparse_sim/) (2026) + Larger projects outside of this repository: diff --git a/cirq_sparse_sim/Demo.ipynb b/cirq_sparse_sim/Demo.ipynb new file mode 100644 index 0000000..1963189 --- /dev/null +++ b/cirq_sparse_sim/Demo.ipynb @@ -0,0 +1,231 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "904fadff-b085-43ed-9a24-f07a892fb238", + "metadata": {}, + "outputs": [], + "source": [ + "### SparseSimulator\n", + "\n", + "This demo compares performance of cirq's built-in dense state simulator (cirq.Simulator) and my SparseSimulator." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "12210309-921d-461f-90eb-2ff6714f2355", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Simulator: n=3, time_elapsed=0.005509\n", + "Simulator: n=4, time_elapsed=0.067351\n", + "Simulator: n=5, time_elapsed=1.071694\n", + "SparseSimulator: n=10, time_elapsed=0.000601\n", + "SparseSimulator: n=100, time_elapsed=0.003026\n", + "SparseSimulator: n=1000, time_elapsed=0.030445\n", + "SparseSimulator: n=10000, time_elapsed=0.429930\n", + "SparseSimulator: n=50000, time_elapsed=4.729522\n" + ] + } + ], + "source": [ + "from collections.abc import Iterator, Sequence\n", + "\n", + "import cirq\n", + "import random\n", + "import time\n", + "\n", + "from sparse_sim import SparseSimulator\n", + "\n", + "\n", + "def maj(a: cirq.Qid, b: cirq.Qid, c: cirq.Qid) -> Iterator[cirq.Operation]:\n", + " yield cirq.CNOT(c, a)\n", + " yield cirq.CNOT(c, b)\n", + " yield cirq.CCNOT(a, b, c)\n", + "\n", + "\n", + "def uma(a: cirq.Qid, b: cirq.Qid, c: cirq.Qid) -> Iterator[cirq.Operation]:\n", + " yield cirq.CCNOT(a, b, c)\n", + " yield cirq.CNOT(c, a)\n", + " yield cirq.CNOT(a, b)\n", + "\n", + "\n", + "class CuccaroAdder(cirq.Gate):\n", + " \"\"\"Adds two n-bit unsigned integers using the Cuccaro ripple-carry circuit.\n", + "\n", + " With little-endian registers, computes\n", + " (a, b, z) -> (a, (a + b) mod 2**n, z XOR floor((a + b) / 2**n)).\n", + "\n", + " Reference: Cuccaro, Draper, Kutin, and Moulton, \"A new quantum\n", + " ripple-carry addition circuit\" (2004), https://arxiv.org/abs/quant-ph/0410184\n", + " \"\"\"\n", + "\n", + " def __init__(self, n: int) -> None:\n", + " if n < 1:\n", + " raise ValueError(\"n must be positive\")\n", + " self.n = n\n", + "\n", + " @property\n", + " def input_sizes(self) -> tuple[int, int, int]:\n", + " # First addend, second addend, carry output.\n", + " return self.n, self.n, 1\n", + "\n", + " def _num_qubits_(self) -> int:\n", + " return 2 * self.n + 1\n", + "\n", + " def _decompose_with_context_(\n", + " self,\n", + " qubits: Sequence[cirq.Qid],\n", + " context: cirq.DecompositionContext,\n", + " ) -> Iterator[cirq.OP_TREE]:\n", + " a = qubits[: self.n]\n", + " b = qubits[self.n : 2 * self.n]\n", + " z = qubits[2 * self.n]\n", + " ancilla = context.qubit_manager.qalloc(1)\n", + " c = ancilla[0]\n", + "\n", + " try:\n", + " yield from maj(c, b[0], a[0])\n", + " for i in range(1, self.n):\n", + " yield from maj(a[i - 1], b[i], a[i])\n", + " yield cirq.CNOT(a[-1], z)\n", + " for i in range(self.n - 1, 0, -1):\n", + " yield from uma(a[i - 1], b[i], a[i])\n", + " yield from uma(c, b[0], a[0])\n", + " finally:\n", + " context.qubit_manager.qfree(ancilla)\n", + "\n", + "\n", + "\n", + "\n", + "rng = random.Random(0)\n", + "simulator = SparseSimulator()\n", + "#simulator = cirq.Simulator()\n", + "\n", + "def run_cuccaro_adder(n, sim): \n", + " a = cirq.LineQubit.range(n)\n", + " b = cirq.LineQubit.range(n, 2 * n)\n", + " z = cirq.LineQubit(2 * n)\n", + " qubits = (*a, *b, z)\n", + " adder = CuccaroAdder(n)\n", + "\n", + " a_value = rng.randint(0, 2**n-1)\n", + " b_value = rng.randint(0, 2**n-1)\n", + " expected_sum = a_value + b_value\n", + "\n", + " circuit = cirq.Circuit(\n", + " cirq.X(qubit)\n", + " for register_value, register in ((a_value, a), (b_value, b))\n", + " for bit, qubit in enumerate(register)\n", + " if (register_value >> bit) & 1\n", + " )\n", + " circuit.append(adder.on(*qubits))\n", + " circuit.append(\n", + " [\n", + " cirq.measure(*a, key=\"a\"),\n", + " cirq.measure(*b, key=\"sum\"),\n", + " cirq.measure(z, key=\"carry_out\"),\n", + " ]\n", + " )\n", + "\n", + " t0 = time.perf_counter()\n", + " result = sim.run(circuit, repetitions=1)\n", + " time_elapsed = time.perf_counter()-t0\n", + " sim_name = sim.__class__.__name__\n", + " print(f\"{sim_name}: {n=}, {time_elapsed=:.6f}\")\n", + " \n", + " measured_a = sum(int(bit) << i for i, bit in enumerate(result.measurements[\"a\"][0]))\n", + " measured_sum = sum(int(bit) << i for i, bit in enumerate(result.measurements[\"sum\"][0]))\n", + " measured_carry = int(result.measurements[\"carry_out\"][0, 0])\n", + " assert measured_a == a_value\n", + " assert measured_sum == expected_sum % 2**n\n", + " assert measured_carry == expected_sum // 2**n\n", + "\n", + "for n in [3,4,5]:\n", + " run_cuccaro_adder(n, cirq.Simulator())\n", + "\n", + "for n in [10, 100, 10**3, 10**4, 5 * 10**4]:\n", + " run_cuccaro_adder(n, SparseSimulator())" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "f13535e7-7bee-431e-a905-d35004c5fc39", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Simulator: n=4, time_elapsed=0.001834\n", + "SparseSimulator: n=4, time_elapsed=0.001021\n", + "Simulator: n=5, time_elapsed=0.002004\n", + "SparseSimulator: n=5, time_elapsed=0.001369\n" + ] + } + ], + "source": [ + "def build_qft_circuit(n, input_value):\n", + " qubits = cirq.LineQubit.range(n)\n", + " circuit = cirq.Circuit(\n", + " cirq.X(qubit)\n", + " for bit, qubit in enumerate(qubits)\n", + " if (input_value >> bit) & 1\n", + " )\n", + " circuit.append(cirq.qft(*qubits))\n", + " return cirq.Circuit(cirq.decompose(circuit)), qubits\n", + "\n", + "\n", + "def benchmark_qft(n, sim):\n", + " input_value = rng.randint(0, 2**n-1)\n", + " circuit, qubits = build_qft_circuit(n, input_value)\n", + " start = time.perf_counter()\n", + " result = simulator.simulate(circuit, qubit_order=qubits)\n", + " time_elapsed = time.perf_counter() - start\n", + "\n", + " sim_name = sim.__class__.__name__\n", + " print(f\"{sim_name}: {n=}, {time_elapsed=:.6f}\")\n", + " return result.final_state_vector\n", + "\n", + "for n in [4,5]:\n", + " benchmark_qft(n, cirq.Simulator())\n", + " benchmark_qft(n, SparseSimulator())\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0c154db6-0088-4f2f-8e82-f5b87e6aa8b0", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From c5bc99bbca9caba7288041c0dfcf0c000ed7bb55 Mon Sep 17 00:00:00 2001 From: Dmytro Fedoriaka Date: Sun, 27 Sep 2026 23:21:25 -0700 Subject: [PATCH 08/11] add QFT test --- .vscode/settings.json | 5 +++ cirq_sparse_sim/Demo.ipynb | 54 ++++++++++++++++-------------- cirq_sparse_sim/sparse_sim_test.py | 34 +++++++++++++++++++ 3 files changed, 67 insertions(+), 26 deletions(-) create mode 100644 .vscode/settings.json diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..1e3a655 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,5 @@ +{ + "chat.tools.terminal.autoApprove": { + ".venv/bin/python": true + } +} \ No newline at end of file diff --git a/cirq_sparse_sim/Demo.ipynb b/cirq_sparse_sim/Demo.ipynb index 1963189..b6e82ba 100644 --- a/cirq_sparse_sim/Demo.ipynb +++ b/cirq_sparse_sim/Demo.ipynb @@ -1,11 +1,9 @@ { "cells": [ { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "id": "904fadff-b085-43ed-9a24-f07a892fb238", "metadata": {}, - "outputs": [], "source": [ "### SparseSimulator\n", "\n", @@ -14,7 +12,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 1, "id": "12210309-921d-461f-90eb-2ff6714f2355", "metadata": {}, "outputs": [ @@ -22,14 +20,14 @@ "name": "stdout", "output_type": "stream", "text": [ - "Simulator: n=3, time_elapsed=0.005509\n", - "Simulator: n=4, time_elapsed=0.067351\n", - "Simulator: n=5, time_elapsed=1.071694\n", - "SparseSimulator: n=10, time_elapsed=0.000601\n", - "SparseSimulator: n=100, time_elapsed=0.003026\n", - "SparseSimulator: n=1000, time_elapsed=0.030445\n", - "SparseSimulator: n=10000, time_elapsed=0.429930\n", - "SparseSimulator: n=50000, time_elapsed=4.729522\n" + "Simulator: n=3, time_elapsed=0.005978\n", + "Simulator: n=4, time_elapsed=0.042375\n", + "Simulator: n=5, time_elapsed=1.041485\n", + "SparseSimulator: n=10, time_elapsed=0.000767\n", + "SparseSimulator: n=100, time_elapsed=0.003140\n", + "SparseSimulator: n=1000, time_elapsed=0.032048\n", + "SparseSimulator: n=10000, time_elapsed=0.427906\n", + "SparseSimulator: n=50000, time_elapsed=4.819097\n" ] } ], @@ -155,7 +153,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 2, "id": "f13535e7-7bee-431e-a905-d35004c5fc39", "metadata": {}, "outputs": [ @@ -163,10 +161,22 @@ "name": "stdout", "output_type": "stream", "text": [ - "Simulator: n=4, time_elapsed=0.001834\n", - "SparseSimulator: n=4, time_elapsed=0.001021\n", - "Simulator: n=5, time_elapsed=0.002004\n", - "SparseSimulator: n=5, time_elapsed=0.001369\n" + "Simulator: n=4, time_elapsed=0.001794\n", + "SparseSimulator: n=4, time_elapsed=0.002759\n", + "Simulator: n=6, time_elapsed=0.002388\n", + "SparseSimulator: n=6, time_elapsed=0.002481\n", + "Simulator: n=8, time_elapsed=0.014807\n", + "SparseSimulator: n=8, time_elapsed=0.006620\n", + "Simulator: n=10, time_elapsed=0.030729\n", + "SparseSimulator: n=10, time_elapsed=0.031616\n", + "Simulator: n=12, time_elapsed=0.139986\n", + "SparseSimulator: n=12, time_elapsed=0.112313\n", + "Simulator: n=14, time_elapsed=0.497365\n", + "SparseSimulator: n=14, time_elapsed=0.507166\n", + "Simulator: n=16, time_elapsed=2.030843\n", + "SparseSimulator: n=16, time_elapsed=2.101261\n", + "Simulator: n=17, time_elapsed=6.624044\n", + "SparseSimulator: n=17, time_elapsed=6.571864\n" ] } ], @@ -193,18 +203,10 @@ " print(f\"{sim_name}: {n=}, {time_elapsed=:.6f}\")\n", " return result.final_state_vector\n", "\n", - "for n in [4,5]:\n", + "for n in [4, 6, 8, 10, 12, 14, 16, 17]:\n", " benchmark_qft(n, cirq.Simulator())\n", " benchmark_qft(n, SparseSimulator())\n" ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "0c154db6-0088-4f2f-8e82-f5b87e6aa8b0", - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { diff --git a/cirq_sparse_sim/sparse_sim_test.py b/cirq_sparse_sim/sparse_sim_test.py index a9809a7..f745467 100644 --- a/cirq_sparse_sim/sparse_sim_test.py +++ b/cirq_sparse_sim/sparse_sim_test.py @@ -161,6 +161,40 @@ def test_single_qubit_unitaries(initial: int, gate: cirq.Gate) -> None: ) +def test_eight_qubit_qft_and_inverse() -> None: + num_qubits = 8 + basis_index = int(np.random.default_rng(0).integers(1 << num_qubits)) + simulator = SparseSimulator() + qubits = simulator.qubit_manager.qalloc(num_qubits) + qft = cirq.QuantumFourierTransformGate(num_qubits) + qft_operation = qft(*qubits) + + transformed = simulator.simulate( + cirq.Circuit(qft_operation), + qubit_order=qubits, + initial_state=basis_index, + ) + + np.testing.assert_allclose( + transformed.final_state_vector, + cirq.unitary(qft)[basis_index], + rtol=0, + atol=1e-9, + ) + + round_trip = simulator.simulate( + cirq.Circuit(qft_operation, cirq.inverse(qft_operation)), + qubit_order=qubits, + ) + + np.testing.assert_allclose( + round_trip.final_state_vector, + np.eye(1 << num_qubits)[0], + rtol=0, + atol=1e-9, + ) + + @pytest.mark.parametrize( "gate", [ From 4603f77272ab52c9ecf347d928f6ebcabaa6a56f Mon Sep 17 00:00:00 2001 From: Dmytro Fedoriaka Date: Mon, 28 Sep 2026 00:21:02 -0700 Subject: [PATCH 09/11] add support for noisy channels. --- cirq_sparse_sim/README.md | 5 +- cirq_sparse_sim/sparse_sim.py | 134 ++++++++++++++++++++++++----- cirq_sparse_sim/sparse_sim_test.py | 121 ++++++++++++++++++++++++-- 3 files changed, 229 insertions(+), 31 deletions(-) diff --git a/cirq_sparse_sim/README.md b/cirq_sparse_sim/README.md index beff4cf..0faf482 100644 --- a/cirq_sparse_sim/README.md +++ b/cirq_sparse_sim/README.md @@ -5,8 +5,9 @@ This is a simple implementation of a Sparse simulator for circuits in Cirq. It is sparse in a sense that it maintains a list of basis states with only non-zero amplitudes. -It supports unitary gates and computational basis measurments, and designed -primarily for simulating arithmetic circuits. +It supports unitary gates, computational basis measurements, and noisy channels, +and is designed primarily for simulating arithmetic circuits. Mixture and Kraus +channels are simulated as stochastic pure-state trajectories. I was not able to find such a simulator in Cirq (maybe I didn't look had enough), so I just wrote this one for my needs. diff --git a/cirq_sparse_sim/sparse_sim.py b/cirq_sparse_sim/sparse_sim.py index 3a2c811..9b72e38 100644 --- a/cirq_sparse_sim/sparse_sim.py +++ b/cirq_sparse_sim/sparse_sim.py @@ -21,7 +21,6 @@ import copy import math import operator -import random from collections.abc import Callable, Iterable, Iterator, Mapping, Sequence from types import NotImplementedType from typing import Any, Self @@ -170,11 +169,11 @@ class _SparseState(cirq.QuantumStateRepresentation, cirq.ClassicalDataStoreReade def __init__( self, - random_source: random.Random, + random_source: np.random.RandomState, qubits: Sequence[cirq.Qid] = (), classical_data: cirq.ClassicalDataStoreReader | None = None, ) -> None: - self.random: random.Random | np.random.RandomState = random_source + self.random = random_source self.axis_by_qubit = {qubit: axis for axis, qubit in enumerate(qubits)} self._next_axis = len(qubits) self._free_axes: list[int] = [] @@ -283,30 +282,104 @@ def _apply_single_qubit_unitary_gate( matrix = cirq.unitary(op_gate, default=None) assert matrix is not None assert matrix.shape == (2, 2) + states, amplitudes = self._apply_matrix(matrix, [qubit_id]) + self._set_normalized_state(states, amplitudes) + + def _apply_matrix( + self, matrix: np.ndarray, qubit_ids: Sequence[int] + ) -> tuple[list[int], list[complex]]: + """Apply a local matrix without normalizing or mutating this state.""" + dimension = 1 << len(qubit_ids) + matrix = np.asarray(matrix) + if matrix.shape != (dimension, dimension): + raise ValueError( + f"Matrix shape {matrix.shape} does not match " + f"{len(qubit_ids)} target qubits" + ) + target_mask = sum(1 << qubit_id for qubit_id in qubit_ids) + nonzero_outputs: dict[int, np.ndarray] = {} merged_amplitudes: dict[int, complex] = {} for state, amplitude in zip(self.basis_states, self.amplitudes, strict=True): - input_bit = (state >> qubit_id) & 1 - for output_bit in (0, 1): - coeff = matrix[output_bit, input_bit] - if coeff == 0: - continue - new_state = ( - state if output_bit == input_bit else state ^ (1 << qubit_id) - ) + input_index = 0 + for qubit_id in qubit_ids: + input_index = (input_index << 1) | ((state >> qubit_id) & 1) + base_state = state & ~target_mask + + outputs = nonzero_outputs.get(input_index) + if outputs is None: + outputs = np.flatnonzero(matrix[:, input_index]) + nonzero_outputs[input_index] = outputs + for output_index in outputs: + new_state = base_state + for position, qubit_id in enumerate(qubit_ids): + output_bit = ( + int(output_index) >> (len(qubit_ids) - position - 1) + ) & 1 + new_state |= output_bit << qubit_id + coeff = matrix[output_index, input_index] merged_amplitudes[new_state] = ( merged_amplitudes.get(new_state, 0.0 + 0.0j) + amplitude * coeff ) + states = list(merged_amplitudes) + amplitudes = list(merged_amplitudes.values()) + return states, amplitudes + + def _set_normalized_state( + self, states: Sequence[int], amplitudes: Sequence[complex] + ) -> None: + """Update sparse state to given states and amplitudes, normalized.""" + norm = math.sqrt(sum(abs(amplitude) ** 2 for amplitude in amplitudes)) + if not math.isfinite(norm) or norm <= 0: + raise ValueError("Selected channel branch has zero or invalid probability") + kept_states = [ - (state, amplitude) - for state, amplitude in merged_amplitudes.items() - if abs(amplitude) >= _AMPLITUDE_EPS + (state, amplitude / norm) + for state, amplitude in zip(states, amplitudes, strict=True) + if abs(amplitude / norm) >= _AMPLITUDE_EPS ] assert len(kept_states) > 0 self.basis_states = [state for state, _ in kept_states] self.amplitudes = _normalize([amplitude for _, amplitude in kept_states]) + def _random_choice(self, probs: list[float]) -> int: + """Sample an index from the given probability distribution.""" + # NumPy rejects empty, negative, non-finite, and non-normalized probabilities. + return int(self.random.choice(len(probs), p=probs)) + + def _apply_mixture( + self, + mixture: Sequence[tuple[float, np.ndarray]], + qubit_ids: Sequence[int], + ) -> None: + """Sample and apply one unitary from a state-independent mixture.""" + probabilities = [float(probability) for probability, _ in mixture] + index = int(self._random_choice(probabilities)) + selected_matrix = mixture[index][1] + + states, amplitudes = self._apply_matrix(np.asarray(selected_matrix), qubit_ids) + self._set_normalized_state(states, amplitudes) + + def _apply_kraus( + self, kraus_operators: Sequence[np.ndarray], qubit_ids: Sequence[int] + ) -> None: + """Sample and normalize one state-dependent Kraus trajectory branch.""" + if not kraus_operators: + raise ValueError("A channel must provide at least one Kraus operator") + + branches = [ + self._apply_matrix(np.asarray(matrix), qubit_ids) + for matrix in kraus_operators + ] + # Each branch probability is the squared norm of K_i |psi>. + probabilities = [ + sum(abs(amplitude) ** 2 for amplitude in amplitudes) + for _, amplitudes in branches + ] + index = self._random_choice(probabilities) + self._set_normalized_state(*branches[index]) + def _measure_qubit(self, target_qubit: int) -> int: """Measure one qubit and collapse sparse state according to Born rule. @@ -431,14 +504,28 @@ def _run_op(self, op: cirq.Operation, context: cirq.DecompositionContext) -> Non return # Prefer explicit unitaries to decompositions introducing unsupported gates. - if ( - op_gate is not None - and len(qubit_ids) == 1 - and cirq.has_unitary(op_gate, allow_decompose=False) - ): + protocol_target = op_gate if op_gate is not None else op + has_unitary = cirq.has_unitary(protocol_target, allow_decompose=False) + if op_gate is not None and len(qubit_ids) == 1 and has_unitary: self._apply_single_qubit_unitary_gate(op_gate, qubit_ids[0]) return + if not has_unitary: + mixture = cirq.mixture(protocol_target, default=None) + if mixture is not None: + self._apply_mixture(mixture, qubit_ids) + return + + kraus_getter = getattr(protocol_target, "_kraus_", None) + if kraus_getter is not None: + kraus_operators = kraus_getter() + if ( + kraus_operators is not None + and kraus_operators is not NotImplemented + ): + self._apply_kraus(tuple(kraus_operators), qubit_ids) + return + num_alloc_before = self.qubit_manager.num_allocated_qubits() decomposed = cirq.decompose_once( op, default=None, context=context, flatten=False @@ -460,8 +547,7 @@ def read_register(self, register: Sequence[cirq.Qid]) -> int: s = self.basis_states[0] try: return sum( - ((s >> self.axis_by_qubit[q]) & 1) << i - for i, q in enumerate(register) + ((s >> self.axis_by_qubit[q]) & 1) << i for i, q in enumerate(register) ) except KeyError as ex: raise ValueError( @@ -638,7 +724,7 @@ class SparseSimulator( def __init__(self, seed: int | None = None) -> None: super().__init__(dtype=np.complex128, seed=seed, split_untangled_states=False) - self.random = random.Random(seed) + self.random = np.random.RandomState(seed) self._state = _SparseState(self.random) self.qubit_manager = SparseSimQubitManager(self) @@ -806,4 +892,6 @@ def _run( pass for key, values in step._classical_data.records.items(): records.setdefault(str(key), []).append(values) - return {key: np.asarray(values, dtype=np.bool_) for key, values in records.items()} + return { + key: np.asarray(values, dtype=np.bool_) for key, values in records.items() + } diff --git a/cirq_sparse_sim/sparse_sim_test.py b/cirq_sparse_sim/sparse_sim_test.py index f745467..613cfa7 100644 --- a/cirq_sparse_sim/sparse_sim_test.py +++ b/cirq_sparse_sim/sparse_sim_test.py @@ -45,6 +45,13 @@ def _prepare_basis( return [cirq.X(q) for i, q in enumerate(qubits) if (value >> i) & 1] +def _random_state_snapshot( + random_state: np.random.RandomState, +) -> tuple[str, bytes, int, int, float]: + algorithm, keys, position, has_gauss, cached_gaussian = random_state.get_state() + return algorithm, keys.tobytes(), position, has_gauss, cached_gaussian + + def _assert_matches_dense( simulator: SparseSimulator, circuit: cirq.Circuit, @@ -1079,8 +1086,6 @@ def _num_qubits_(self) -> int: "gate", [ _UnsupportedGate(), - cirq.depolarize(0.1), - cirq.amplitude_damp(0.1), ], ids=repr, ) @@ -1092,6 +1097,110 @@ def test_unsupported_operations_raise_value_error(gate: cirq.Gate) -> None: simulator.run(cirq.Circuit(gate(*qubits))) +@pytest.mark.parametrize( + "channel, initial, expected_one_probability", + [ + (cirq.bit_flip(0.25), 0, 0.25), + ( + cirq.asymmetric_depolarize(p_x=0.1, p_y=0.2, p_z=0.3), + 0, + 0.3, + ), + (cirq.depolarize(0.3), 0, 0.2), + (cirq.amplitude_damp(0.3), 1, 0.7), + (cirq.generalized_amplitude_damp(p=0.7, gamma=0.4), 1, 0.72), + ], + ids=repr, +) +def test_common_channels_follow_expected_basis_state_statistics( + channel: cirq.Gate, initial: int, expected_one_probability: float +) -> None: + simulator = SparseSimulator(seed=1234) + qubits = simulator.qubit_manager.qalloc(1) + result = simulator.run( + cirq.Circuit( + _prepare_basis(qubits, initial), + channel(*qubits), + cirq.measure(*qubits, key="m"), + ), + repetitions=400, + ) + + assert np.mean(result.measurements["m"]) == pytest.approx( + expected_one_probability, abs=0.08 + ) + + +def test_phase_flip_trajectory_has_expected_x_statistics() -> None: + probability = 0.25 + simulator = SparseSimulator(seed=1234) + qubits = simulator.qubit_manager.qalloc(1) + result = simulator.run( + cirq.Circuit( + cirq.H(*qubits), + cirq.phase_flip(probability)(*qubits), + cirq.H(*qubits), + cirq.measure(*qubits, key="m"), + ), + repetitions=400, + ) + + assert np.mean(result.measurements["m"]) == pytest.approx(probability, abs=0.08) + + +def test_phase_damping_trajectory_has_expected_x_statistics() -> None: + gamma = 0.36 + simulator = SparseSimulator(seed=1234) + qubits = simulator.qubit_manager.qalloc(1) + result = simulator.run( + cirq.Circuit( + cirq.H(*qubits), + cirq.phase_damp(gamma)(*qubits), + cirq.H(*qubits), + cirq.measure(*qubits, key="m"), + ), + repetitions=400, + ) + + expected_one_probability = (1 - math.sqrt(1 - gamma)) / 2 + assert np.mean(result.measurements["m"]) == pytest.approx( + expected_one_probability, abs=0.08 + ) + + +def test_multiqubit_depolarizing_channel_uses_cirq_qubit_order() -> None: + simulator = SparseSimulator(seed=1234) + qubits = simulator.qubit_manager.qalloc(2) + result = simulator.run( + cirq.Circuit( + cirq.depolarize(1.0, n_qubits=2)(*qubits), + cirq.measure(*qubits, key="m"), + ), + repetitions=600, + ) + + outcomes = np.array( + [cirq.big_endian_bits_to_int(bits) for bits in result.measurements["m"]] + ) + frequencies = np.bincount(outcomes, minlength=4) / len(outcomes) + np.testing.assert_allclose(frequencies, [3 / 15, 4 / 15, 4 / 15, 4 / 15], atol=0.08) + + +def test_channel_trajectories_are_seeded() -> None: + circuit = cirq.Circuit( + cirq.H(cirq.LineQubit(0)), + cirq.amplitude_damp(0.4)(cirq.LineQubit(0)), + cirq.measure(cirq.LineQubit(0), key="m"), + ) + + results = [ + SparseSimulator(seed=123).run(circuit, repetitions=100).measurements["m"] + for _ in range(2) + ] + + np.testing.assert_array_equal(results[0], results[1]) + + def test_unsupported_classical_condition_raises_value_error() -> None: simulator = SparseSimulator() qubits = simulator.qubit_manager.qalloc(1) @@ -1262,7 +1371,7 @@ def test_zero_repetitions_preserve_state_results_and_random_stream() -> None: simulator = SparseSimulator(seed=5) qubits = simulator.qubit_manager.qalloc(2) simulator.run(cirq.Circuit(cirq.X(qubits[0]), cirq.measure(*qubits, key="old"))) - random_state = simulator.random.getstate() + random_state = _random_state_snapshot(simulator.random) circuit = cirq.Circuit( cirq.H(qubits[0]), cirq.CircuitOperation(cirq.FrozenCircuit(cirq.measure(*qubits, key="new"))), @@ -1273,7 +1382,7 @@ def test_zero_repetitions_preserve_state_results_and_random_stream() -> None: assert result.measurements["new"].shape == (0, 2) assert simulator.basis_states == [1] assert simulator.measurement_results == {"old": 1} - assert simulator.random.getstate() == random_state + assert _random_state_snapshot(simulator.random) == random_state @pytest.mark.parametrize("repetitions", [1, 3]) @@ -1414,7 +1523,7 @@ def test_step_sampling_is_seeded_and_does_not_collapse_state(repetitions: int) - ) step = list(simulator.simulate_moment_steps(circuit))[-1] before = step.state_vector() - random_state = simulator.random.getstate() + random_state = _random_state_snapshot(simulator.random) samples = step.sample([qubits[2], qubits[1], qubits[0]], repetitions, seed=42) @@ -1426,7 +1535,7 @@ def test_step_sampling_is_seeded_and_does_not_collapse_state(repetitions: int) - ) np.testing.assert_allclose(step.state_vector(), before) assert step.measurements == {} - assert simulator.random.getstate() == random_state + assert _random_state_snapshot(simulator.random) == random_state assert step.sample([], repetitions, seed=42).shape == (repetitions, 0) with pytest.raises(ValueError, match="repetitions"): step.sample(qubits, repetitions=-1) From 51d96b11df5844aa58992d2557ee8e0f38b40eae Mon Sep 17 00:00:00 2001 From: Dmytro Fedoriaka Date: Mon, 28 Sep 2026 00:22:09 -0700 Subject: [PATCH 10/11] remove conversion --- cirq_sparse_sim/sparse_sim.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cirq_sparse_sim/sparse_sim.py b/cirq_sparse_sim/sparse_sim.py index 9b72e38..fd88ff4 100644 --- a/cirq_sparse_sim/sparse_sim.py +++ b/cirq_sparse_sim/sparse_sim.py @@ -355,7 +355,7 @@ def _apply_mixture( ) -> None: """Sample and apply one unitary from a state-independent mixture.""" probabilities = [float(probability) for probability, _ in mixture] - index = int(self._random_choice(probabilities)) + index = self._random_choice(probabilities) selected_matrix = mixture[index][1] states, amplitudes = self._apply_matrix(np.asarray(selected_matrix), qubit_ids) From f84aadb4ddde96327dca3610801ef1628e159442 Mon Sep 17 00:00:00 2001 From: Dmytro Fedoriaka Date: Mon, 28 Sep 2026 00:25:05 -0700 Subject: [PATCH 11/11] add .vscode to gitignore --- .gitignore | 3 ++- .vscode/settings.json | 5 ----- 2 files changed, 2 insertions(+), 6 deletions(-) delete mode 100644 .vscode/settings.json diff --git a/.gitignore b/.gitignore index 1f12c0a..afa0c38 100644 --- a/.gitignore +++ b/.gitignore @@ -6,4 +6,5 @@ local .idea .pytest_cache obj -.venv \ No newline at end of file +.venv +.vscode \ No newline at end of file diff --git a/.vscode/settings.json b/.vscode/settings.json deleted file mode 100644 index 1e3a655..0000000 --- a/.vscode/settings.json +++ /dev/null @@ -1,5 +0,0 @@ -{ - "chat.tools.terminal.autoApprove": { - ".venv/bin/python": true - } -} \ No newline at end of file