Running Simulations

Execute circuits: sf.run(), sf.simulate(), 4 simulation methods, sf.State, noise models, and RunResult.

Two Entry Points

sf.run() returns a RunResult (counts + state + metadata). sf.simulate() returns sf.State directly.

import superfermion as sf

qc = sf.Circuit(2).h(0).cnot(0, 1)

# Run: counts + state + metadata
result = sf.run(qc, device="cpu", shots=1024)
print(result.counts)   # {'00': 512, '11': 512}
print(result.state)    # sf.State

# Simulate: state only (shots=0 implied)
state = sf.simulate(qc, device="cpu")
print(state.numpy())   # [0.707+0j, 0, 0, 0.707+0j]

Devices

# Local simulation
sf.run(qc, device="cpu")    # Rust CPU — Rayon multithreading
sf.run(qc, device="gpu")    # Rust GPU — CUDA (requires sf-gpu)

# Cloud QPUs
from superfermion.devices.ibm import IBMDevice
from superfermion.devices.ionq import IonQDevice

ibm = IBMDevice(token="...")
sf.run(qc, device=ibm("ibm_brisbane"), shots=8192)

ionq = IonQDevice(api_key="...")
sf.run(qc, device=ionq("aria-1"), shots=1000)

Simulation Methods

Change method= — the API doesn't change.

Statevector (default)

Exact simulation. Up to ~25 qubits CPU, ~30 GPU.

result = sf.run(qc, method="statevector", shots=1024)
state = sf.simulate(qc, method="statevector")

MPS (Tensor Network)

For low-entanglement circuits. Scales to 200+ qubits.

result = sf.run(qc, method="mps", bond_dim=64, shots=1024)

Higher bond dimension = more accuracy, more memory. Start at 64.

MPS also executes dynamic (mid-circuit) circuits with shots>0: each shot replays the circuit while collapsing the tensor network at every mid measure, with bond_dim capping the bond (default 64; reported in result.metadata["bond_dim"]). See Mid-Circuit Measurement & Feed-Forward.

MPS shots=0 densifies the tensor network into a full statevector. For more than 26 qubits this raises a MemoryError with guidance — pass shots>0 to sample directly from the tensor network instead, or densify=False to skip densification and receive only the MPS State handle (result.state.numpy() densifies once on demand). MPS sampling and statevector output use the same qubit order as the statevector backend.

Stabilizer

Clifford-only circuits. Polynomial time, scales to ~4000 qubits.

clifford = sf.Circuit(100).h(0)
for i in range(99):
    clifford = clifford.cnot(i, i + 1)

result = sf.run(clifford, method="stabilizer", shots=1024)

Non-Clifford gates (T, RZ, etc.) raise an error. Use statevector for those.

State operations the stabilizer backend does not support — .numpy(), .partial_trace(), .fidelity(), .entropy(), .purity(), .probabilities(), .grad() — raise superfermion.MethodError (a RuntimeError subclass).

Density Matrix

Noisy simulation with Kraus channels. Up to ~12 qubits.

noise = sf.NoiseModel().add_depolarizing(0.001).add_amplitude_damping(0.002)
result = sf.run(qc, method="density_matrix", noise_model=noise, shots=0)

Key and layout conventions:

counts / probabilities keys are little-endian — qubit q is bit q of int(key, 2) — identical to statevector, mps and the default sampling path (x(2) on 3 qubits reports "100").
result.metadata["density_matrix"] is a 2-D array in big-endian order (qubit 0 = most significant): for x(0) on 2 qubits the weight sits at rho[2, 2].
sf.State.numpy() returns the density matrix as a flat little-endian array (length 4^n); reshape to (2**n, 2**n) to read the matrix (the bit-reversed image of the metadata rho).

Method comparison

Method	Max Qubits	Exact	Supports Gradients	Supports Noise
`statevector`	~25	Yes	Yes	No
`mps`	200+	Approximate	No	No
`stabilizer`	~4000	Yes (Clifford)	No	No
`density_matrix`	~12	Yes	No	Yes

`sf.State` — the Rust quantum state

sf.State wraps a Rust quantum state object. All methods dispatch to compiled Rust.

state = sf.simulate(qc)

# Basic properties
print(state.n_qubits)       # 2
print(state.shape)          # (4,)
print(state.method)         # "statevector"
print(state.device)         # "cpu"

# Extract data
sv = state.numpy()          # complex ndarray
samples = state.sample(1000)  # list[int] — measurement outcomes

# Derived quantities
print(state.entropy())          # Von Neumann entropy (0 = pure)
print(state.purity())           # 1.0 = pure, <1 = mixed
print(state.probabilities())    # Outcome probability distribution

# Expectation values
obs = [([3, 3], 1.0, 0.0)]    # ZZ observable
energy = state.expectation(obs)
print(energy)                  # 1.0

# Gradients — on a parameterized ansatz
ansatz = sf.Circuit(2).ry(sf.param("t"), 0).cnot(0, 1)
state_p = sf.simulate(ansatz, params={"t": 0.5})
dag = ansatz.to_ir()  # keep parameters symbolic — .bind() would erase them
params = {"t": 0.5}
grads = state_p.grad(obs, dag, params)
print(grads)                   # {'t': 0.0}

# Operations on states
print(state.fidelity(state))         # 1.0 (self-fidelity)
rdm = state.partial_trace([0])       # Reduced density matrix
qfim = state_p.qfim(dag, params)     # Quantum Fisher Information Matrix
print(qfim)

# Create from numpy
import numpy as np
vec = np.array([1, 0, 0, 0], dtype=complex)
state2 = sf.State.from_numpy(vec, n_qubits=2)
print(state2.numpy())

`RunResult` — execution result

obs = [([3, 3], 1.0, 0.0)]    # ZZ observable
result = sf.run(qc, device="cpu", shots=8192)

# Primary outputs
print(result.counts)           # {'00': 4096, '11': 4096}
print(result.state)            # sf.State
print(result.shots)            # 8192

# Derived
print(result.probabilities)        # {'00': 0.5, '11': 0.5} (all backends)
print(result.get_probabilities())    # {'00': 0.5, '11': 0.5}
print(result.statevector)      # ndarray
print(result.circuit)          # original circuit
print(result.metadata)         # execution metadata

# Methods
energy = result.expectation(obs)       # expectation from state
print(energy)
result.plot()                          # histogram plot (saves to file)

Noise Models

noise = (sf.NoiseModel()
    .add_depolarizing(0.001)           # symmetric depolarizing
    .add_amplitude_damping(0.002)      # T1 relaxation
    .add_phase_damping(0.001)          # pure dephasing
    .add_bit_flip(0.0005)              # bit flip
    .add_two_qubit_depolarizing(0.001) # correlated two-qubit noise
)

qc = sf.Circuit(1).x(0)
result = sf.run(qc, method="density_matrix", noise_model=noise, shots=0)
print(f"purity = {result.metadata['purity']:.4f}")

Auto-Bind Parameters

ansatz = sf.Circuit(2).ry(sf.param("t"), 0).cnot(0, 1)

# Both forms are equivalent:
state = sf.simulate(ansatz, params={"t": 0.5})
state = sf.simulate(ansatz.bind({"t": 0.5}))

Error Handling

from superfermion import MethodError

try:
    # T gate on a stabilizer circuit
    sf.run(sf.Circuit(1).t(0), method="stabilizer")
except RuntimeError as e:
    # Non-Clifford gates cannot be simulated with the stabilizer method.
    # MethodError is a RuntimeError subclass, so `except RuntimeError`
    # also catches stabilizer methods like .numpy() / .partial_trace().
    print(f"Method not applicable: {e}")

Bit Order & Conventions

Statevectors are little-endian: qubit q lives at bit position q of the amplitude index, so state.numpy()[3] is the amplitude of |q0=1, q1=1⟩.
Bitstring keys are q0-last: in counts and probabilities, "001" means qubit 0 is |1⟩ and qubits 1, 2 are |0⟩. This matches format(index, f"0{n}b") for the little-endian amplitude index.
Observable labels are q0-leftmost: in SparsePauliOp.from_string("ZI"), the first character acts on qubit 0 — "ZI" is Z⊗I (Z on qubit 0).
counts, probabilities, and state.numpy() all agree with each other and with the expectation-value paths (sf.expval, state.expectation, adjoint gradients).

On this page