Python Bindings

prism-q ships Python bindings built with PyO3. They are a thin wrapper over the Rust crate: the compiled extension is prism_q._prism_q and the pure-Python prism_q package re-exports it. Simulation runs in Rust with the GIL released, so the wrapper adds no per-gate overhead.

Wheels are abi3 for Python 3.11 and newer, so one wheel per platform covers every supported interpreter.

Install

pip install prism-q

NumPy is the only runtime dependency. Building from a source checkout needs maturin:

pip install maturin
maturin develop --manifest-path bindings/python/Cargo.toml

The bindings enable the parallel feature by default. The gpu feature is optional and off in the published wheels (see GPU backends); the distributed backend is not reachable from Python.

Quick start

from prism_q import CircuitBuilder, simulate

circuit = CircuitBuilder(2, 2).h(0).cx(0, 1).measure_all().build()
counts = simulate(circuit).seed(42).shots(1000).counts()
print(counts)          # {'00': 507, '11': 493}

Building circuits

CircuitBuilder is a fluent API. Every gate method returns the builder, and build() produces the Circuit that simulation consumes.

from prism_q import CircuitBuilder

circuit = (
    CircuitBuilder(3, 3)
    .h(0)
    .cx(0, 1)
    .rz(0.5, 2)
    .cphase(0.25, 1, 2)
    .measure_all()
    .build()
)

pauli_rotation(theta, factors) appends exp(-i * theta * P / 2) for the Pauli string P given as (qubit, axis) factors with axis one of "X", "Y", "Z"; identity factors are omitted. A weight-1 string lowers to rx, ry, or rz and a two-qubit ZZ string to rzz, so fusion and Clifford recognition keep firing on them. Circuit.add_pauli_rotation is the imperative spelling.

builder.pauli_rotation(0.4, [(0, "X"), (1, "Y"), (3, "Z")]).param(0)

cu and mcu take a 2x2 matrix as nested Python sequences of complex numbers. Out-of-range qubits raise PrismError at build time rather than at simulation time.

Three other routes produce a Circuit:

from prism_q import Circuit, circuits, parse_qasm

manual = Circuit(2, 2)                    # imperative, add_gate / add_measure / add_reset
ghz = circuits.ghz(10)                    # pre-built corpus
parsed = parse_qasm(qasm_source)          # OpenQASM 3.0, with 2.0 accepted

The circuits submodule mirrors the Rust builders documented in Circuit Builders: qft, ghz, w_state, random, hardware_efficient_ansatz, clifford_heavy, clifford_random_pairs, qaoa, single_qubit_rotation, clifford_t, quantum_volume, cz_chain, phase_estimation, independent_bell_pairs, independent_random_blocks, and local_clifford_blocks. Seeded builders default to seed 42.

Running a simulation

simulate(circuit) returns a Simulation you configure with .seed(), .backend(), and .noise(), then finish with a terminal method. The default seed is 42.

from prism_q import BackendKind, simulate

sim = simulate(circuit).seed(7).backend(BackendKind.statevector())
outcome = sim.run()

shots() and sample_counts() average trajectories on any backend holding a per-shot pure state. The three rows marked "density matrix only" read the exact mixed state instead, so they need .backend(BackendKind.density_matrix()); auto dispatch never selects it. There the mixture is evolved once and every terminal reads that one evolution, so the probabilities are seed independent and the observables carry no sampling error. Circuits with mid-circuit measurement or classical conditioning are rejected on that route, since the mixture holds every measurement branch at once.

Terminals that cannot honor a model raise PrismError naming the reason, rather than silently ignoring it. state_vector() always uses the statevector backend and density_matrix_expectation_values() always uses the density-matrix backend, both regardless of .backend(...).

ShotsResult and CountsResult both expose counts(), returning a dict keyed by bitstring.

Distributions too wide to write down

A circuit whose qubits fall into independent groups is answered per group, and RunOutcome keeps it that way: the dense vector is built only when probabilities is read. Fifteen independent Bell pairs span 30 qubits, whose dense form is 8 GB, and the blocks are 15 arrays of four entries.

outcome = simulate(circuits.independent_bell_pairs(15)).seed(42).run()

outcome.num_basis_states            # 2 ** 30, and nothing was materialized
for qubits, probs in outcome.probabilities_factored():
    print(qubits, probs)            # [0, 1] [0.5 0. 0. 0.5], ...

Each block is (qubits, probs) with qubits ascending, and probs indexed by those qubits packed in that order with qubits[0] in the least significant bit. The probability of a basis state is the product of one entry per block, which is what probabilities computes.

probabilities_factored() returns None when the run produced a dense distribution, which is the common case: the decomposed route needs the widest group several qubits narrower than the register, so two Bell pairs stay dense. num_basis_states is None when the backend exposed no distribution at all.

Result metadata

RunOutcome, ShotsResult, and CountsResult each carry a metadata object describing how the result was produced.

result = simulate(circuit).seed(42).run()
print(result.metadata.backend)               # 'Statevector'
print(result.metadata.engine)                # None unless samplers share the backend
print(result.metadata.is_exact)              # True
print(result.metadata.fidelity_lower_bound)  # None when exact
print(result.metadata.placement)             # 'host' or 'device'

is_exact is False when the engine that ran can discard state weight or estimate by sampling. It marks the route, not the run: an MPS whose bond dimension the circuit never fills reports is_exact == False with fidelity_lower_bound == 1.0, so the flag answers whether the answer could have been approximated and the bound answers whether it was.

Automatic dispatch sends a circuit past the statevector cap to a bounded-bond MPS, which is the only route those circuits have. That is taken by default and the result says so. .require_exact() rejects it instead, raising PrismError naming the engine it would have used.

simulate(big_circuit).seed(42).require_exact().marginals()  # raises

Starting from a state other than |0...0>

.initial_state(amplitudes) replaces the default all-zero start. It takes any sequence of complex numbers, a complex128 NumPy array included, indexed with qubit 0 in the least significant bit.

import math
import numpy as np
from prism_q import CircuitBuilder, simulate

theta = math.pi / 8
start = np.array([math.cos(theta), math.sin(theta)], dtype=np.complex128)
circuit = CircuitBuilder(1).h(0).build()
probs = simulate(circuit).initial_state(start).seed(42).run().probabilities

The vector must have 2 ** num_qubits entries and unit norm. A wrong length, a non-finite entry, or a norm off unity raises PrismError; an unnormalized vector is rejected rather than rescaled, so a mistake surfaces instead of becoming a silent factor on every amplitude.

A start state also narrows the route. Auto dispatch reads circuit structure, and its shortcuts (tableau, product state, subsystem decomposition, Pauli propagation) are only valid from |0...0>: a Clifford circuit produces a stabilizer state only when its input is one. So auto() resolves to the statevector, the GPU and distributed statevectors and density_matrix() are the only other backends that accept one, and every other choice raises PrismError naming itself. run(), shots(), sample_counts(), marginals(), expectation_values(), and state_vector() carry it; expectation_gradient() and density_matrix_expectation_values() reject it, as do shots() and sample_counts() with a noise model attached, since trajectory replay reinitializes a pure state per shot. To evolve a start state under noise, read the exact mixture with run(), marginals(), or expectation_values() on density_matrix().

Selecting a backend

BackendKind.auto() is the default and picks a backend from circuit structure. Pass an explicit one to override it.

The density-matrix backend stores 4^n amplitudes, so its qubit ceiling is about half the statevector cap; exceeding it raises PrismError naming the cap. The distributed statevector backend has no Python constructor: MPI_Init ownership between the interpreter, mpi4py, and the extension is unsettled.

GPU backends

The GPU constructors take a GpuContext, an opaque handle to one CUDA device and its compiled kernels. Build it once and reuse it: construction compiles the kernel module, and passing the same handle to several simulations shares that work.

from prism_q import BackendKind, GpuContext, circuits, simulate

context = GpuContext(0)
outcome = simulate(circuits.qft(16)).backend(BackendKind.auto_gpu(context)).seed(42).run()
print(outcome.probabilities)

GpuContext(device_id) is where a missing or unusable device is reported, and it raises PrismError rather than falling back. Past construction, routing is soft by design and matches the Rust API: statevector_gpu runs circuits below the crossover (PRISM_GPU_MIN_QUBITS, default 14) on the host, auto_gpu routes each block independently, and a block whose device allocation fails degrades to the host rather than erroring. A run that produces host results is therefore normal, not a failure signal.

stabilizer_gpu sets its crossover at 100000 qubits (PRISM_STABILIZER_GPU_MIN_QUBITS), so it runs on the host tableau unless that override is lowered. The device tableau is correct; the default stays high until benchmarks justify lowering it.

density_matrix_gpu(context) holds the exact mixed state on the device and is the one device kind with no crossover and no host fallback: auto_gpu never selects it, every noisy terminal that density_matrix() serves answers from the device buffer, and a width whose 4^n buffer does not fit in free device memory raises PrismError before anything is allocated (13 qubits on an 11 GiB card).

The published wheels are built without CUDA, because two of the three wheel targets have no CUDA toolkit and macOS has no CUDA at all. In those wheels the constructors still exist and GpuContext(...) raises PrismError naming the missing build feature, so code written against the GPU API fails with a message rather than an AttributeError. Two predicates separate the cases:

GpuContext.is_supported()   # was this build compiled with CUDA support
GpuContext.is_available()   # ... and is a usable device present

To get a build with CUDA support, install the CUDA toolkit (12.x or newer) and build from a checkout:

maturin develop --manifest-path bindings/python/Cargo.toml --features gpu

On Windows that build links the toolkit's NVRTC library (nvrtc64_120_0.dll for CUDA 12.x) from the toolkit bin directory, which Python does not search. The package adds it on import when CUDA_PATH is set, which the toolkit installer does; without it the import fails with DLL load failed while importing _prism_q.

Distributed backend

The distributed statevector shards the dense state across MPI ranks. It is reachable from Python through a DistributedContext, which attaches to an MPI that is already running. This extension never calls MPI_Init or MPI_Finalize: mpi4py does both, at its own import and at interpreter exit, and a handle whose refcount drop finalized MPI would make every later MPI call in the process erroneous.

from mpi4py import MPI  # starts MPI; import before touching the context
from prism_q import BackendKind, DistributedContext, circuits, simulate

context = DistributedContext()
outcome = (
    simulate(circuits.qft(24))
    .backend(BackendKind.statevector_distributed(context))
    .seed(42)
    .run()
)
print(context.rank, context.size, outcome.probabilities[:4])

Run it with mpiexec -n 4 python script.py. The rank count must be a power of two, and every rank needs enough local qubits after the shard bits are taken (PRISM_DIST_MIN_LOCAL_QUBITS).

The contract is SPMD, and it is enforced rather than assumed. Four ranks are four interpreters running the same source, and every collective inside the backend is entered by all of them, so a script that branches before the call

if comm.rank == 0:
    result = simulate(circuit).backend(...).run()   # deadlocks

hangs the job: rank 0 blocks in a collective the others never enter. Every rank calls run with the same circuit and seed; only what you do with the returned value may branch on rank. Two cross-checks turn the common violations into errors instead of hangs, each costing one collective at run entry: ranks disagreeing about the register shape, the seed, or the tuning knobs are rejected, and so are ranks handed different circuits.

Two surfaces are deliberately loud rather than convenient.

A world of one rank raises. MPI-2 and later make a singleton MPI_Init succeed, so a script launched without mpiexec would otherwise get a correct answer from one rank at single-host speed with no signal that nothing was distributed. Pass allow_single_rank=True when that is what you meant.

Constructing the context without MPI running raises rather than starting MPI, so a forgotten from mpi4py import MPI is reported at the point it happened.

The thread level mpi4py negotiated must be at least MPI_THREAD_FUNNELED: Rayon workers run beside MPI in the same process, and every MPI call stays on the thread that constructed the context. mpi4py requests MPI_THREAD_MULTIPLE by default; a script that lowers it (mpi4py.rc.thread_level = 'single' or mpi4py.rc.threads = False) gets an error from the constructor rather than a run that mixes threads with a single-threaded MPI.

The published wheels have no MPI support: mpi-sys runs bindgen and needs a system MPI at build time, and the extension has to link the same MPI implementation and ABI as mpi4py and as the launcher. Mixing two implementations in one process corrupts rather than failing loudly, which is why this is a from-source path:

maturin develop --manifest-path bindings/python/Cargo.toml --features distributed-mpi

DistributedContext.is_supported() says whether a build has MPI support, the same way GpuContext.is_supported() does for CUDA.

Sub-communicators, per-rank device placement, and any distributed path other than the statevector backend are out of scope: the context is the world communicator and nothing else.

Noise

Build a NoiseModel from a circuit, then attach it. A model is sized to the circuit it was built from, and using it with a different circuit raises PrismError.

from prism_q import NoiseChannel, NoiseModel, simulate

model = NoiseModel.uniform_depolarizing(circuit, 0.01)
counts = simulate(circuit).seed(42).noise(model).sample_counts(4000).counts()

NoiseModel.uniform_depolarizing(circuit, p) and NoiseModel.with_amplitude_damping(circuit, gamma) cover the common cases. For per-instruction control, start from NoiseModel.empty(circuit) and attach events:

model = NoiseModel.empty(circuit)
model.add_event(0, NoiseChannel.amplitude_damping(0.05), [0])
model.add_event(1, NoiseChannel.two_qubit_depolarizing(0.02), [0, 1])
model.with_readout_error(0.01, 0.01)
model.validate()

Channels are pauli(px, py, pz), depolarizing(p), amplitude_damping(gamma), phase_damping(gamma), thermal_relaxation(t1, t2, gate_time), two_qubit_depolarizing(p), and custom(kraus) for an explicit list of 2x2 Kraus operators. validate() checks probabilities and Kraus completeness; is_pauli_only() reports whether the model holds only single-qubit Pauli channels and no readout error. A model carrying readout error or two_qubit_depolarizing answers False there and still runs on the stabilizer samplers, which apply readout to the measurement record and sample the pair channel as one joint draw over its 15 branches.

Expectation values

An observable is a list of (qubit, axis) factors with axis one of "X", "Y", "Z". Identity factors are omitted, so [(0, "Z"), (2, "X")] means Z0 ⊗ I1 ⊗ X2. Both expectation terminals take a list of observables and return one float each.

observables = [[(0, "Z")], [(0, "Z"), (1, "Z")]]
values = simulate(circuit).seed(42).expectation_values(observables)

expectation_values requires a unitary circuit and gives ⟨ψ|P|ψ⟩. density_matrix_expectation_values evolves the density matrix through the circuit and any attached noise model and gives exact Tr(rho P), with measurements read off the final mixed state without collapse. It is the zero-variance analogue of averaging over trajectories:

model = NoiseModel.empty(circuit)
model.add_event(0, NoiseChannel.amplitude_damping(0.3), [0])
exact = simulate(circuit).seed(42).noise(model).density_matrix_expectation_values(
    [[(0, "Z")]]
)

Parameter sweeps

A variational loop rebinds angles while the gate sequence stays fixed. Parameters names the slots those angles land in, and PreparedCircuit holds the circuit across bindings so fusion and backend selection are settled once rather than per point.

from prism_q import CircuitBuilder, PreparedCircuit

builder = CircuitBuilder(4)
for q in range(4):
    builder.ry(0.1, q).param(q)
for q in range(3):
    builder.cx(q, q + 1)

prepared = PreparedCircuit(builder.build(), builder.parameters())
for values in [[0.1, 0.2, 0.3, 0.4], [0.5, 0.6, 0.7, 0.8]]:
    outcome = prepared.run(values, seed=42)

run(values, seed=42) returns the same RunOutcome as simulate(...).run(). bind(values) returns the bound Circuit instead, for handing to simulate(...) with different options or to any other consumer of a circuit.

Automatic dispatch reads the template, so build it at angles representative of the sweep. A template whose rotations are all zero reads as Clifford and settles on a backend that then rejects the bound circuit. Pass an explicit backend as the third argument to decide it yourself:

prepared = PreparedCircuit(circuit, params, BackendKind.statevector())

reuses_fusion_plan reports whether the recorded fused structure is being replayed. It is a performance fact, not an error: results agree either way.

Build a Parameters three ways. builder.parameters() returns the set recorded by param(slot); Parameters.all_rotations(circuit) gives every bindable gate its own slot in circuit order; Parameters(n) plus link(instruction, slot) declares the slots up front. Several gates may share a slot, in which case binding writes one angle to each.

A wrong-length value vector, a non-finite angle, and a link pointing at a gate that carries no angle all raise PrismError.

Gradients

expectation_gradient computes ⟨H⟩ and its exact gradient with respect to the bound parameters by the adjoint method, at a cost independent of the parameter count. Mark parameters with param(slot) while building, then pass parameter_links() through.

import numpy as np
from prism_q import CircuitBuilder, simulate

builder = CircuitBuilder(2)
builder.ry(0.3, 0).param(0)
builder.cx(0, 1)
builder.rz(0.7, 1).param(1)
circuit, links = builder.build(), builder.parameter_links()

hamiltonian = [(1.0, [(0, "Z")]), (0.5, [(0, "Z"), (1, "Z")])]
value, gradient = simulate(circuit).seed(42).expectation_gradient(hamiltonian, links)

A Hamiltonian term is (coefficient, observable). Several gates may share a slot, in which case their gradients accumulate. param() rejects anything but a differentiable gate (rx, ry, rz, rzz, p, pauli_rotation), and the circuit must be unitary.

Quantum error correction

QecProgram exposes the native QEC IR: reset, measure, detector, observable_include, postselect, and noise, with QecBasis, QecNoise, and RecordRef as the supporting types. run() returns a QecResult carrying detector, observable, and measurement arrays as NumPy bool_ matrices, plus logical_error_rates() and survivor_rate(). Programs can also be parsed from text with QecProgram.from_text. See the Noise and QEC guide for the model itself.

detector_error_model() derives the program's DetectorErrorModel for decoding: probabilities() (float64), detector_matrix(), and observable_matrix() (bool, detectors or observables by mechanisms) feed check-matrix decoder constructors directly, detector_coords() carries the per-detector coordinates, decompose_graphlike() returns the form matching decoders need (at most two detectors per mechanism), and to_text() writes the common detector error model text format for file-based decoders:

dem = qp.detector_error_model()
H, p, L = dem.detector_matrix(), dem.probabilities(), dem.observable_matrix()
with open("memory_d3.dem", "w") as f:
    f.write(dem.to_text())

Decoder runs the built-in union-find decoder over a graphlike model, so a memory experiment's logical error rate needs no external decoder. decode takes the (shots, num_detectors) bool detector array and returns the (shots, num_observables) predicted observable flips:

decoder = prism_q.Decoder(dem.decompose_graphlike())
res = qp.run()
predicted = decoder.decode(res.detectors)
failures = (predicted[:, 0] != res.observables[:, 0]).sum()

Errors and typing

Every failure surfaces as prism_q.PrismError, carrying the message from the Rust error. Backend limits, unsupported operations, and invalid arguments all raise it rather than panicking:

import prism_q

try:
    simulate(huge).density_matrix_expectation_values([[(0, "Z")]])
except prism_q.PrismError as exc:
    print(exc)   # backend `density_matrix` is incompatible: circuit has ... qubits

The package ships type stubs (prism_q/_prism_q.pyi) and a py.typed marker, so mypy and Pyright see the full surface.