Adds tools for appending randomized measurement bases and processing renyi entropy from bitstring

cirq-core/cirq/experiments/__init__.py

@@ -71,3 +71,8 @@
     parallel_two_qubit_xeb,
     run_rb_and_xeb,
 )
+
+from cirq.experiments.measure_in_random_bases import (
+    append_randomized_measurements,
+    RandomizedMeasurements,
+)

cirq-core/cirq/experiments/measure_in_random_bases.py

@@ -0,0 +1,143 @@
+# Copyright 2024 The Cirq Developers
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     https://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from collections.abc import Mapping, Sequence
+from typing import Any, Literal
+
+import cirq
+import numpy as np
+import numpy.typing as npt
+
+
+class RandomizedMeasurements:
+    def __init__(
+        self,
+        num_qubits: int,
+        num_unitaries: int,
+        subsystem: Sequence[str | int] | None = None,
+        qubit_mapping: Mapping[int, str | int] | None = None,
+        rng: np.random.Generator = np.random.default_rng(),
+    ):
+        """Class structure for performing and analyzing a general randomized measurement protocol.
+        For more details on the randomized measurement toolbox see https://arxiv.org/abs/2203.11374
+        Args:
+            num_qubits: Number of qubits in the circuit
+            num_unitaries: Number of random pre-measurement unitaries
+            subsystem: The specific subsystem measured in random basis
+            qubit_mapping: The mapping between the measurement bitstring index and qubit specifier
+            rng: Random generator
+        """
+        self.num_qubits = num_qubits
+        self.num_unitaries = num_unitaries
+        self.subsystem = subsystem
+
+        self.qubit_mapping = (
+            qubit_mapping if qubit_mapping else {i: i for i in range(self.num_qubits)}
+        )
+
+        self.rng = rng
+
+        self.pre_measurement_unitaries_list = self._generate_unitaries_list()
+
+    def _generate_unitaries_list(self) -> npt.NDArray[Any]:
+        """Generates a list of pre-measurement unitaries."""
+
+        pauli_strings = self.rng.choice(["X", "Y", "Z"], size=(self.num_unitaries, self.num_qubits))
+
+        if self.subsystem is not None:
+            for i in range(pauli_strings.shape[1]):
+                if self.qubit_mapping[i] not in self.subsystem:
+                    pauli_strings[:, i] = np.array(["Z"] * self.num_unitaries)
+
+        return pauli_strings
+
+    def unitaries_to_moment(
+        self, unitaries: Sequence[Literal["X", "Y", "Z"]], qubits: Sequence[Any]
+    ) -> 'cirq.Moment':
+        """Outputs the cirq moment associated with the pre-measurement rotations.
+        Args:
+            unitaries: List of pre-measurement unitaries
+            qubits: List of qubits
+
+        Returns: The cirq moment associated with the pre-measurement rotations
+        """
+        op_list: list[cirq.Operation] = []
+        for idx, pauli in enumerate(unitaries):
+            op_list.append(_pauli_basis_rotation(pauli).on(qubits[idx]))
+
+        return cirq.Moment.from_ops(*op_list)
+
+
+def _pauli_basis_rotation(basis: Literal["X", "Y", "Z"]) -> 'cirq.Gate':
+    """Given a measurement basis returns the associated rotation.
+    Args:
+        basis: Measurement basis
+    Returns: The cirq gate for associated with measurement basis
+    """
+    if basis == "X":
+        return cirq.Ry(rads=-np.pi / 2)
+    elif basis == "Y":
+        return cirq.Rx(rads=np.pi / 2)
+    elif basis == "Z":
+        return cirq.I
+
+
+def append_randomized_measurements(
+    circuit: 'cirq.AbstractCircuit',
+    randomized_measurements_generator: RandomizedMeasurements | None = None,
+    *,
+    subsystem: tuple[int] | None = None,
+    qubits: Sequence | None = None,
+    num_unitaries: int | None = None,
+    rng: np.random.Generator = np.random.default_rng(),
+) -> Sequence['cirq.Circuit']:
+    """Given an input circuit returns a list of circuits with the pre-measurement unitaries.
+    If no arguments are specified, it will default to computing the entropy of the entire
+    circuit.
+
+    Args:
+        circuit: The input circuit
+        randomized_measurements_generator: RandomizedMeasurements class to use for
+        generating random measurements.
+        subsystem: The specific subsystem measured in random basis.
+        qubits: A sequence of qubits to measure in random basis.
+        num_unitaries: The number of random pre-measurement unitaries to append.
+        rng: Random number genergate
+    Returns:
+        List of circuits with pre-measurement unitaries and measurements added
+    """
+    qubits = qubits or list(circuit.all_qubits())
+
+    if randomized_measurements_generator is None:
+        randomized_measurements_generator = RandomizedMeasurements(
+            len(qubits),
+            num_unitaries if num_unitaries else len(qubits),
+            subsystem=subsystem,
+            rng=rng,
+        )
+
+    circuit_list = []
+
+    for unitaries in randomized_measurements_generator.pre_measurement_unitaries_list:
+        pre_measurement_moment = randomized_measurements_generator.unitaries_to_moment(
+            unitaries, qubits
+        )
+
+        temp_circuit = cirq.Circuit.from_moments(
+            *circuit.moments, pre_measurement_moment, cirq.measure_each(*qubits)
+        )
+
+        circuit_list.append(temp_circuit)
+
+    return circuit_list

cirq-core/cirq/experiments/measure_in_random_bases_test.py

@@ -0,0 +1,119 @@
+# Copyright 2024 The Cirq Developers
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     https://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import cirq
+import cirq.experiments.measure_in_random_bases as mrb
+
+
+def test_append_randomized_measurements_generates_n_circuits():
+    # Create a 4-qubit circuit
+    q0, q1, q2, q3 = cirq.LineQubit.range(4)
+    circuit = cirq.Circuit([cirq.H(q0), cirq.CNOT(q0, q1), cirq.CNOT(q1, q2), cirq.CNOT(q2, q3)])
+
+    # Append randomized measurements
+    circuits = mrb.append_randomized_measurements(circuit)
+
+    assert len(circuits) == 4  # num of qubits
+
+
+def test_append_randomized_measurements_appends_two_moments_to_end_of_circuit():
+    # Create a 4-qubit circuit
+    q0, q1, q2, q3 = cirq.LineQubit.range(4)
+    circuit = cirq.Circuit([cirq.H(q0), cirq.CNOT(q0, q1), cirq.CNOT(q1, q2), cirq.CNOT(q2, q3)])
+
+    num_moments_pre = len(circuit.moments)
+
+    # Append randomized measurements
+    circuits = mrb.append_randomized_measurements(circuit)
+
+    for circuit in circuits:
+        num_moments_post = len(circuit.moments)
+        assert num_moments_post == num_moments_pre + 2  # 1 random gate and 1 measurement gate
+
+
+def test_append_randomized_measurements_generates_n_circuits_when_passed_subystem_arg():
+    # Create a 4-qubit circuit
+    q0, q1, q2, q3 = cirq.LineQubit.range(4)
+    circuit = cirq.Circuit([cirq.H(q0), cirq.CNOT(q0, q1), cirq.CNOT(q1, q2), cirq.CNOT(q2, q3)])
+
+    # Append randomized measurements to subsystem
+    circuits = mrb.append_randomized_measurements(circuit, subsystem=(0, 1))
+
+    assert len(circuits) == 4
+
+
+def test_append_randomized_measurements_leaves_qubits_not_in_specified_subsystem_unchanged():
+    # Create a 4-qubit circuit
+    q0, q1, q2, q3 = cirq.LineQubit.range(4)
+    circuit = cirq.Circuit([cirq.H(q0), cirq.CNOT(q0, q1), cirq.CNOT(q1, q2), cirq.CNOT(q2, q3)])
+
+    # Append randomized measurements to subsystem
+    circuits = mrb.append_randomized_measurements(circuit, subsystem=(0, 1))
+
+    for circuit in circuits:
+        # assert latter subsystems were not changed.
+        assert circuit.operation_at(q2, 4) == cirq.I(q2)
+        assert circuit.operation_at(q3, 4) == cirq.I(q3)
+
+
+def test_append_randomized_measurements_generates_circuits_only_for_passed_qubit_mapping():
+    # Create a 4-qubit circuit
+    q0, q1, q2, q3 = cirq.LineQubit.range(4)
+    circuit = cirq.Circuit([cirq.H(q0), cirq.CNOT(q0, q1), cirq.CNOT(q1, q2), cirq.CNOT(q2, q3)])
+
+    # Append randomized measurements to subsystem
+    circuits = mrb.append_randomized_measurements(circuit, qubits=[q0, q1])
+
+    assert len(circuits) == 2
+
+
+def test_append_randomized_measurements_appends_two_moments_to_specified_qubits():
+    # Create a 4-qubit circuit
+    q0, q1, q2, q3 = cirq.LineQubit.range(4)
+    circuit = cirq.Circuit([cirq.H(q0), cirq.CNOT(q0, q1), cirq.CNOT(q1, q2), cirq.CNOT(q2, q3)])
+    num_moments_pre = len(circuit.moments)
+
+    # Append randomized measurements to subsystem
+    circuits = mrb.append_randomized_measurements(circuit, qubits=[q0, q1])
+
+    for circuit in circuits:
+        num_moments_post = len(circuit.moments)
+        assert num_moments_post == num_moments_pre + 2  # 1 random gate and 1 measurement gate
+
+
+def test_append_randomized_measurements_with_num_unitaries_generates_k_circuits():
+    # Create a 4-qubit circuit
+    q0, q1, q2, q3 = cirq.LineQubit.range(4)
+    circuit = cirq.Circuit([cirq.H(q0), cirq.CNOT(q0, q1), cirq.CNOT(q1, q2), cirq.CNOT(q2, q3)])
+    num_unitaries = 3
+
+    # Append randomized measurements to subsystem
+    circuits = mrb.append_randomized_measurements(circuit, num_unitaries=num_unitaries)
+
+    assert len(circuits) == num_unitaries
+
+
+def test_append_randomized_measurements_with_num_unitaries_appends_two_moments_on_each_circuit():
+    # Create a 4-qubit circuit
+    q0, q1, q2, q3 = cirq.LineQubit.range(4)
+    circuit = cirq.Circuit([cirq.H(q0), cirq.CNOT(q0, q1), cirq.CNOT(q1, q2), cirq.CNOT(q2, q3)])
+    num_moments_pre = len(circuit.moments)
+    num_unitaries = 3
+
+    # Append randomized measurements to subsystem
+    circuits = mrb.append_randomized_measurements(circuit, num_unitaries=num_unitaries)
+
+    for circuit in circuits:
+        num_moments_post = len(circuit.moments)
+        assert num_moments_post == num_moments_pre + 2

cirq-core/cirq/qis/__init__.py

@@ -60,3 +60,4 @@
     average_error,
     decoherence_pauli_error,
 )
+from cirq.qis.process_renyi_entropy_from_bitstrings import process_entropy_from_bitstrings

cirq-core/cirq/qis/process_renyi_entropy_from_bitstrings.py

@@ -0,0 +1,116 @@
+# Copyright 2024 The Cirq Developers
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     https://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from concurrent.futures import ThreadPoolExecutor
+from collections.abc import Sequence
+from typing import Any, Optional
+
+import numpy as np
+import numpy.typing as npt
+
+
+def _get_hamming_distance(
+    bitstring_1: npt.NDArray[np.int8], bitstring_2: npt.NDArray[np.int8]
+) -> int:
+    """Calculates the Hamming distance between two bitstrings.
+    Args:
+        bitstring_1: Bitstring 1
+        bitstring_2: Bitstring 2
+    Returns: The Hamming distance
+    """
+    return (bitstring_1 ^ bitstring_2).sum().item()
+
+
+def _bitstrings_to_probs(
+    bitstrings: npt.NDArray[np.int8],
+) -> tuple[npt.NDArray[np.int8], npt.NDArray[Any]]:
+    """Given a list of bitstrings from different measurements returns a probability distribution.
+    Args:
+        bitstrings: The bitstring
+    Returns:
+        A tuple of bitstrings and their corresponding probabilities.
+    """
+
+    num_shots = bitstrings.shape[0]
+    unique_bitstrings, counts = np.unique(bitstrings, return_counts=True, axis=0)
+    probs = counts / num_shots
+
+    return (unique_bitstrings, probs)
+
+
+def _bitstring_format_helper(
+    measured_bitstrings: npt.NDArray[np.int8], subsystem: Sequence[int] | None = None
+) -> npt.NDArray[np.int8]:
+    """Formats the bitstring for analysis based on the selected subsystem.
+    Args:
+        measured_bitstrings: Measured bitstring
+        subsystem: Subsystem of interest
+    Returns: The bitstring string for the subsystem
+    """
+    if subsystem is None:
+        return measured_bitstrings
+
+    return measured_bitstrings[:, :, subsystem]
+
+
+def _compute_bitstring_purity(bitstrings: npt.NDArray[np.int8]) -> float:
+    """Computes the purity of a bitstring.
+    Args:
+        bitstrings: The bitstrings measured using the same unitary operators
+    Returns: The purity of the bitstring
+    """
+
+    bitstrings, probs = _bitstrings_to_probs(bitstrings)
+    purity = 0
+    for idx, s in enumerate(bitstrings):
+        p = probs[idx]
+        for j, s_prime in enumerate(bitstrings):
+            p_prime = bitstrings[j]
+            purity += (-2.0) ** float(-_get_hamming_distance(s, s_prime)) * p * p_prime
+
+    return purity * 2 ** (bitstrings.shape[-1])
+
+
+def process_entropy_from_bitstrings(
+    measured_bitstrings: npt.NDArray[np.int8],
+    subsystem: tuple[int] | None = None,
+    pool: Optional[ThreadPoolExecutor] = None,
+) -> float:
+    """Compute the renyi entropy of an array of bitstrings.
+    Args:
+        measured_bitstrings: List of numpy arrays.
+        subsystem: Subsystem of interest
+        pool: ThreadPoolExecutor used to paralelleize the computation.
+
+    Returns:
+        A float indicating the computed entropy.
+    """
+    bitstrings = _bitstring_format_helper(measured_bitstrings, subsystem)
+    num_shots = bitstrings.shape[1]
+    num_qubits = bitstrings.shape[-1]
+
+    if num_shots == 1:
+        return 0
+
+    if pool is not None:
+        with pool as executor:
+            purities = list(executor.map(_compute_bitstring_purity, list(bitstrings)))
+        purity = np.mean(purities)
+
+    else:
+        purity = np.mean([_compute_bitstring_purity(bitstring) for bitstring in bitstrings])
+
+    purity_unbiased = purity * num_shots / (num_shots - 1) - (2**num_qubits) / (num_shots - 1)
+
+    return purity_unbiased