"""CuttingPrimitive: PrimitiveAlgorithm for QPD circuit cutting.
build() uses streaming experiment generation (O(n_cuts × n_gates) peak memory)
and QWC-grouped sub-blocks (one block per QWC group per experiment × subcircuit
instead of one per observable term).
run() post-processes engine SamplingResults using bitmask parity — one result
feeds all terms in a QWC group simultaneously.
Set max_experiments to run a random subset of the 6^n_cuts experiments for an
approximate (but unbiased) result at reduced cost.
"""
from math import ceil
from typing import Optional
import numpy as np
import qarpx as qx
from qarp._types import Shots
from qarp.algorithms._primitives.primitive_algorithm import PrimitiveAlgorithm
from qarp.algorithms._primitives.target import Target
from qarp.blocks import AnyBlock
from qarp.cutting import EAPartitioning, QPDDecomposition, Reconstructer, check_cut_budget
from qarp.cutting._qpd_decomposition import _post_process_shot
from qarp.errors import CapabilityError
from qarp.operators import QubitOperator
[docs]
class CuttingPrimitive(PrimitiveAlgorithm):
"""PrimitiveAlgorithm that evaluates expectation values via QPD circuit cutting.
Works with any OpenQARP engine without modification.
"""
supported_targets = frozenset({Target.EXPECTATION_VALUE})
# QPD experiment sampling is inherently stochastic — no ∞-shot limit.
supports_exact = False
# The QPD reconstruction is linear in each fragment's expectation value taken
# separately — all a per-occurrence shift needs (products across fragments are fine).
gradient_kind = "expectation"
def __init__(
self,
ket: Optional[AnyBlock] = None,
operator: Optional[QubitOperator] = None,
max_subcircuit_qubits: Optional[int] = None,
n_shots: int = 5000,
force_max_number_cuts: bool = False,
experiment_fraction: Optional[float] = None,
rng_seed: Optional[int] = None,
):
"""
Args:
ket: Quantum state block to measure.
operator: Observable (``qarp.operators.QubitOperator``).
max_subcircuit_qubits: Capacity cap per subcircuit (qubit count
each partition must fit into); defaults to 256.
n_shots: Shots per sub-experiment.
force_max_number_cuts: Skip the config.max_number_of_cuts guard.
experiment_fraction: Fraction of the 6^n_cuts experiments to run,
as a value in ``(0, 1]``. ``None`` (default) runs all
experiments. The sampling strategy is selected automatically:
``"top_k"`` for parametric cut gates (RZZ etc.) and
``"uniform"`` for non-parametric gates (CX). Inspect
``prim._qpd.sampling_strategy`` after ``build()`` to see which
was chosen.
rng_seed: Seed for the experiment-selection RNG (``"uniform"``
mode only).
"""
if n_shots is Shots.EXACT:
# supports_exact = False; EXACT would also slip past the
# `n_shots is not None` assert in build().
raise CapabilityError(
"CuttingPrimitive has no ∞-shot limit (QPD experiment sampling "
"is inherently stochastic); pass a finite n_shots."
)
super().__init__(
ket=ket,
operator=operator,
n_shots=n_shots,
target=Target.EXPECTATION_VALUE,
)
self.max_subcircuit_qubits = max_subcircuit_qubits or 256
self.force_max_number_cuts = force_max_number_cuts
self.experiment_fraction = experiment_fraction
self.rng_seed = rng_seed
[docs]
def circuits(self) -> list:
"""Partition the circuit and build all measurement sub-blocks.
The cut reads the ket's gate stream, so the ket is built here, one
of two primitives (with ``PauliShadow``) whose circuits cannot be
derived from an unbuilt input. Uses streaming experiment generation and QWC grouping so
that the circuits number n_experiments × n_subcircuits ×
n_qwc_groups (instead of × n_obs_terms).
"""
if self.ket is None:
raise RuntimeError("ket must be set before calling build().")
if self.operator is None:
raise RuntimeError("operator must be set before calling build().")
# 1. Flatten the ket block, strip barriers and measurements
built_ket = self.ket.build() if not self.ket.is_built else self.ket
raw_commands = list(built_ket.flatten())
n_qubits = built_ket.n_qubits
rec = Reconstructer(raw_commands, n_qubits)
commands, _ = rec.remove_measure_gates(
[c for c in raw_commands if c.gate != qx.GateType.Barrier]
)
# 2. Partition into subcircuits
max_n = self.max_subcircuit_qubits
n_subcircuits = ceil(n_qubits / max_n)
max_sizes = [max_n] * n_subcircuits
self._cutter = EAPartitioning(commands, n_qubits, max_sizes, verbose=False)
self._cutter_result = self._cutter.cut()
check_cut_budget(self._cutter_result.n_cuts, force=self.force_max_number_cuts)
# 3. QPD decomposition — stream experiments + QWC block building
# CuttingPrimitive always carries an int n_shots (defaulted in __init__),
# though the base attribute is typed Optional[Union[int, Shots]].
assert isinstance(self.n_shots, int) # CuttingPrimitive rejects Shots.EXACT in __init__
self._qpd = QPDDecomposition(
self._cutter_result,
observable=self.operator,
verbose=False,
n_shots=self.n_shots,
experiment_fraction=self.experiment_fraction,
rng_seed=self.rng_seed,
force_max_number_cuts=self.force_max_number_cuts,
)
# Split observable and compute QWC groups (does NOT materialise experiments)
self._qpd._split_observable()
# Build blocks streaming: yields one experiment at a time, builds QWC groups
circuits, self._metadata = self._qpd.build_experiment_blocks_streaming()
return circuits
[docs]
def run(self, results: list) -> float:
"""Reconstruct the expectation value from engine sampling results.
Args:
results: list[qx.SamplingResult] in the same order as sub_blocks.
Returns:
Reconstructed expectation value (float).
"""
n_exps = len(self._qpd.coefficients)
# shots_data[str(i_sub)][i_exp][term_idx] = ev
shots_data: dict = {
str(i): [{} for _ in range(n_exps)] for i in range(self._qpd.n_subcircuits)
}
for result, (
k_exp,
i_sub,
_grp_idx,
term_indices,
qpd_cbits,
obs_cbits,
local_masks,
) in zip(results, self._metadata, strict=True):
local_n = len(obs_cbits)
if len(result.cbit_history) == 0:
if qpd_cbits or obs_cbits:
# The experiment recorded measurements but the backend
# returned no per-shot register — treating that as "all
# contributions are 1" would be silently wrong.
raise RuntimeError(
"CuttingPrimitive: experiment expects classical bits "
f"(qpd={list(qpd_cbits)}, obs={list(obs_cbits)}) but the "
"SamplingResult has an empty cbit_history — the engine's "
"simulator does not fill the per-shot classical register."
)
# No measurements recorded — all contributions are 1
for term_idx in term_indices:
shots_data[str(i_sub)][k_exp][term_idx] = 1.0
else:
accum = {term_idx: 0.0 for term_idx in term_indices}
for shot_cbits in result.cbit_history:
qpd_vals = [int(shot_cbits[i]) for i in qpd_cbits]
obs_bits_int = sum(int(shot_cbits[obs_cbits[k]]) << k for k in range(local_n))
qpd_factor = _post_process_shot(qpd_vals)
for term_idx, mask in zip(term_indices, local_masks, strict=True):
pauli_factor = (-1) ** bin(obs_bits_int & mask).count("1")
accum[term_idx] += qpd_factor * pauli_factor
n = result.n_shots
for term_idx in term_indices:
shots_data[str(i_sub)][k_exp][term_idx] = accum[term_idx] / n
exp_vals = self._qpd._reconstruct_expectation_value_from_shots(shots_data)
return float(
np.real(
np.dot(exp_vals, self._qpd.sub_observables_coeffs) + self._qpd.observable.constant
)
)