Draft support for reconstructing a probability distribution

circuit_knitting/cutting/cutting_decomposition.py

@@ -115,18 +115,11 @@ def cut_gates(
     Returns:
         A copy of the input circuit with the specified gates replaced with :class:`.TwoQubitQPDGate`\ s
         and a list of :class:`.QPDBasis` instances -- one for each decomposed gate.
-
-    Raises:
-        ValueError: The input circuit should contain no classical bits or registers.
     """
-    if len(circuit.cregs) != 0 or circuit.num_clbits != 0:
-        raise ValueError(
-            "Circuits input to cut_gates should contain no classical registers or bits."
-        )
-    # Replace specified gates with TwoQubitQPDGates
     if not inplace:
         circuit = circuit.copy()
 
+    # Replace specified gates with TwoQubitQPDGates
     bases = []
     for gate_id in gate_ids:
         gate = circuit.data[gate_id]

circuit_knitting/cutting/cutting_experiments.py

@@ -82,14 +82,21 @@ def generate_cutting_experiments(
             to the same cut.
         ValueError: :class:`SingleQubitQPDGate` instances are not allowed in unseparated circuits.
     """
-    if isinstance(circuits, QuantumCircuit) and not isinstance(observables, PauliList):
-        raise ValueError(
-            "If the input circuits is a QuantumCircuit, the observables must be a PauliList."
-        )
-    if isinstance(circuits, dict) and not isinstance(observables, dict):
-        raise ValueError(
-            "If the input circuits are contained in a dictionary keyed by partition labels, the input observables must also be represented by such a dictionary."
-        )
+    if observables is None:
+        # FIXME: ensure there's at least one measurement in at least one of the subsystems.
+        # And it's kind of weird, because any subcircuit without measurements does not need to be run!!
+        pass
+    else:
+        if isinstance(circuits, QuantumCircuit) and not isinstance(
+            observables, PauliList
+        ):
+            raise ValueError(
+                "If the input circuits is a QuantumCircuit, the observables must be a PauliList."
+            )
+        if isinstance(circuits, dict) and not isinstance(observables, dict):
+            raise ValueError(
+                "If the input circuits are contained in a dictionary keyed by partition labels, the input observables must also be represented by such a dictionary."
+            )
     if not num_samples >= 1:
         raise ValueError("num_samples must be at least 1.")
 
@@ -99,13 +106,16 @@ def generate_cutting_experiments(
     if isinstance(circuits, QuantumCircuit):
         is_separated = False
         subcircuit_dict: dict[Hashable, QuantumCircuit] = {"A": circuits}
-        subobservables_by_subsystem = decompose_observables(
-            observables, "A" * len(observables[0])
-        )
-        subsystem_observables = {
-            label: ObservableCollection(subobservables)
-            for label, subobservables in subobservables_by_subsystem.items()
-        }
+        if observables is None:
+            subsystem_observables = None
+        else:
+            subobservables_by_subsystem = decompose_observables(
+                observables, "A" * len(observables[0])
+            )
+            subsystem_observables = {
+                label: ObservableCollection(subobservables)
+                for label, subobservables in subobservables_by_subsystem.items()
+            }
         # Gather the unique bases from the circuit
         bases, qpd_gate_ids = _get_bases(circuits)
         subcirc_qpd_gate_ids: dict[Hashable, list[list[int]]] = {"A": qpd_gate_ids}
@@ -120,10 +130,12 @@ def generate_cutting_experiments(
         bases = _get_bases_by_partition(subcircuit_dict, subcirc_qpd_gate_ids)
 
         # Create the commuting observable groups
-        subsystem_observables = {
-            label: ObservableCollection(so) for label, so in observables.items()
-        }
-
+        if observables is None:
+            subsystem_observables = None
+        else:
+            subsystem_observables = {
+                label: ObservableCollection(so) for label, so in observables.items()
+            }
     # Sample the joint quasiprobability decomposition
     random_samples = generate_qpd_weights(bases, num_samples=num_samples)
 
@@ -144,17 +156,23 @@ def generate_cutting_experiments(
         sampled_coeff = (redundancy / num_samples) * (kappa * np.sign(actual_coeff))
         coefficients.append((sampled_coeff, weight_type))
         map_ids_tmp = map_ids
-        for label, so in subsystem_observables.items():
-            subcircuit = subcircuit_dict[label]
+        for label, subcircuit in subcircuit_dict.items():
             if is_separated:
                 map_ids_tmp = tuple(map_ids[j] for j in subcirc_map_ids[label])
-            for j, cog in enumerate(so.groups):
-                new_qc = _append_measurement_register(subcircuit, cog)
-                decompose_qpd_instructions(
-                    new_qc, subcirc_qpd_gate_ids[label], map_ids_tmp, inplace=True
+            if subsystem_observables is None:
+                new_qc = decompose_qpd_instructions(
+                    subcircuit, subcirc_qpd_gate_ids[label], map_ids_tmp, inplace=False
                 )
-                _append_measurement_circuit(new_qc, cog, inplace=True)
                 subexperiments_dict[label].append(new_qc)
+            else:
+                so = subsystem_observables[label]
+                for j, cog in enumerate(so.groups):
+                    new_qc = _append_measurement_register(subcircuit, cog)
+                    decompose_qpd_instructions(
+                        new_qc, subcirc_qpd_gate_ids[label], map_ids_tmp, inplace=True
+                    )
+                    _append_measurement_circuit(new_qc, cog, inplace=True)
+                    subexperiments_dict[label].append(new_qc)
 
     # Remove initial and final resets from the subexperiments.  This will
     # enable the `Move` operation to work on backends that don't support

circuit_knitting/cutting/cutting_reconstruction.py

@@ -24,11 +24,31 @@
 
 from ..utils.observable_grouping import CommutingObservableGroup, ObservableCollection
 from ..utils.bitwise import bit_count
+from ..utils.iteration import strict_zip
 from .cutting_decomposition import decompose_observables
 from .cutting_experiments import _get_pauli_indices
 from .qpd import WeightType
 
 
+def reconstruct_distribution(
+    results: SamplerResult,
+    original_circuit_num_clbits: int,
+    coefficients: Sequence[tuple[float, WeightType]],
+):
+    """Reconstruct probability distribution."""
+    num_meas_bits = original_circuit_num_clbits
+    quasi_dists_out: dict[str | int, float] = {}
+    for quasi_dist, (coeff, _) in strict_zip(results.quasi_dists, coefficients):
+        for outcome, weight in quasi_dist.items():
+            meas_outcomes = outcome & ((1 << num_meas_bits) - 1)
+            qpd_outcomes = outcome >> num_meas_bits
+            qpd_factor = 1 - 2 * (bit_count(qpd_outcomes) & 1)
+            quasi_dists_out[meas_outcomes] = (
+                quasi_dists_out.get(meas_outcomes, 0.0) + coeff * qpd_factor * weight
+            )
+    return quasi_dists_out
+
+
 def reconstruct_expectation_values(
     results: (
         SamplerResult

test/cutting/test_cutting_decomposition.py

@@ -266,15 +266,6 @@ def test_cut_gates(self):
             qc.cx(0, 1)
             qpd_qc, _ = cut_gates(qc, [0])
             self.assertEqual(qpd_qc, compare_qc)
-        with self.subTest("classical bit on input"):
-            qc = QuantumCircuit(2, 1)
-            qc.cx(0, 1)
-            with pytest.raises(ValueError) as e_info:
-                cut_gates(qc, [0])
-            assert (
-                e_info.value.args[0]
-                == "Circuits input to cut_gates should contain no classical registers or bits."
-            )
 
     def test_unused_qubits(self):
         """Issue #218"""

test/cutting/test_cutting_roundtrip.py

@@ -15,8 +15,8 @@
 import logging
 
 import numpy as np
-from qiskit import QuantumCircuit
-from qiskit.circuit.library import UnitaryGate
+from qiskit.circuit import QuantumCircuit, ClassicalRegister
+from qiskit.circuit.library import EfficientSU2, UnitaryGate
 from qiskit.circuit.library.standard_gates import (
     RXXGate,
     RYYGate,
@@ -47,9 +47,11 @@
 from circuit_knitting.utils.simulation import ExactSampler
 from circuit_knitting.cutting import (
     partition_problem,
+    cut_gates,
     generate_cutting_experiments,
     reconstruct_expectation_values,
 )
+from circuit_knitting.cutting.cutting_reconstruction import reconstruct_distribution
 from circuit_knitting.cutting.instructions import Move
 
 logger = logging.getLogger(__name__)
@@ -224,3 +226,42 @@ def test_sampler_with_identity_subobservable(sampler, is_exact_sampler):
 
         # Ensure both methods yielded equivalent expectation values
         assert np.allclose(exact_expvals, reconstructed_expvals, atol=1e-8)
+
+
+@pytest.fixture(
+    params=[
+        (8, 4, 1, 0.4),
+    ]
+)
+def example_sampler_circuit(request):
+    num_qubits, num_measurements, reps, gate_param = request.param
+
+    circuit0 = EfficientSU2(
+        num_qubits=num_qubits, reps=reps, entanglement="circular"
+    ).decompose()
+    circuit0.assign_parameters([gate_param] * len(circuit0.parameters), inplace=True)
+
+    circuit0.add_register(ClassicalRegister(num_measurements))
+    for i in range(num_measurements):
+        circuit0.measure(i, i)
+
+    cut_indices = [
+        i
+        for i, instruction in enumerate(circuit0.data)
+        if {circuit0.find_bit(q)[0] for q in instruction.qubits} == {0, num_qubits - 1}
+    ]
+    circuit1, bases = cut_gates(circuit0, cut_indices)
+
+    return circuit0, circuit1, bases
+
+
+def test_cutting_exact_distribution_reconstruction(example_sampler_circuit):
+    circuit0, circuit1, bases = example_sampler_circuit
+    subexperiments, coefficients = generate_cutting_experiments(circuit1, None, np.inf)
+    subexperiment_results = ExactSampler().run(subexperiments).result()
+    reconstructed = reconstruct_distribution(
+        subexperiment_results, circuit1.num_clbits, coefficients
+    )
+    exact = ExactSampler().run(circuit0).result().quasi_dists[0]
+    for k in range(2**circuit1.num_clbits):
+        assert reconstructed.get(k, 0.0) == pytest.approx(exact.get(k, 0.0))