Mermin refactor

src/qibocal/protocols/__init__.py

@@ -72,6 +72,7 @@
     chsh_pulses,
     correct_virtual_z_phases,
     correct_virtual_z_phases_signal,
+    mermin,
     optimize_two_qubit_gate,
 )
 from .two_qubit_state_tomography import two_qubit_state_tomography
@@ -149,5 +150,6 @@
     "standard_rb_2q",
     "standard_rb_2q_inter",
     "optimize_two_qubit_gate",
+    "mermin",
     "ramsey_zz",
 ]

src/qibocal/protocols/two_qubit_interaction/__init__.py

@@ -1,5 +1,6 @@
 from .chevron import chevron, chevron_signal
 from .chsh import chsh_circuits, chsh_pulses
+from .mermin import mermin
 from .optimize import optimize_two_qubit_gate
 from .virtual_z_phases import correct_virtual_z_phases
 from .virtual_z_phases_signal import correct_virtual_z_phases_signal

src/qibocal/protocols/two_qubit_interaction/mermin/__init__.py

@@ -0,0 +1 @@
+from .protocol import mermin

src/qibocal/protocols/two_qubit_interaction/mermin/protocol.py

@@ -0,0 +1,259 @@
+from dataclasses import dataclass, field
+from typing import Optional
+
+import numpy as np
+import numpy.typing as npt
+import plotly.graph_objects as go
+from qibolab import ExecutionParameters
+from qibolab.platform import Platform
+from qibolab.qubits import QubitId
+
+from qibocal.auto.operation import Data, Parameters, Results, Routine
+
+from ...readout_mitigation_matrix import readout_mitigation_matrix
+from ...utils import STRING_TYPE, calculate_frequencies
+from .pulses import create_mermin_sequences
+from .utils import (
+    compute_mermin,
+    get_mermin_coefficients,
+    get_mermin_polynomial,
+    get_readout_basis,
+)
+
+
+@dataclass
+class MerminParameters(Parameters):
+    """Mermin experiment input parameters."""
+
+    ntheta: int
+    """Number of angles probed linearly between 0 and 2 pi."""
+    native: Optional[bool] = False
+    """If True a circuit will be created using only GPI2 and CZ gates."""
+    apply_error_mitigation: Optional[bool] = False
+    """Error mitigation model"""
+
+
+MerminType = np.dtype(
+    [
+        ("theta", float),
+        ("basis", STRING_TYPE),
+        ("state", STRING_TYPE),
+        ("frequency", int),
+    ]
+)
+
+
+@dataclass
+class MerminData(Data):
+    """Mermin Data structure."""
+
+    thetas: list
+    """Angles probed."""
+    data: dict[list[QubitId], npt.NDArray[MerminType]] = field(default_factory=dict)
+    """Raw data acquired."""
+    mitigation_matrix: dict[list[QubitId], npt.NDArray[np.float64]] = field(
+        default_factory=dict
+    )
+    """Mitigation matrix computed using the readout_mitigation_matrix protocol."""
+
+    @property
+    def targets(self):
+        return list(self.data)
+
+
+@dataclass
+class MerminResults(Results):
+    """Mermin Results class."""
+
+    mermin: dict[tuple[QubitId, ...], npt.NDArray[np.float64]] = field(
+        default_factory=dict
+    )
+    """Raw Mermin value."""
+
+    mermin_mitigated: dict[tuple[QubitId, ...], npt.NDArray[np.float64]] = field(
+        default_factory=dict
+    )
+    """Mitigated Mermin value."""
+
+
+def _acquisition(
+    params: MerminParameters,
+    platform: Platform,
+    targets: list[list[QubitId]],
+) -> MerminData:
+    """Data acquisition for Mermin protocol using pulse sequences."""
+
+    thetas = np.linspace(0, 2 * np.pi, params.ntheta)
+    data = MerminData(thetas=thetas.tolist())
+    if params.apply_error_mitigation:
+        mitigation_data, _ = readout_mitigation_matrix.acquisition(
+            readout_mitigation_matrix.parameters_type.load(
+                dict(pulses=True, nshots=params.nshots)
+            ),
+            platform,
+            targets,
+        )
+
+        mitigation_results, _ = readout_mitigation_matrix.fit(mitigation_data)
+        data.mitigation_matrix = mitigation_results.readout_mitigation_matrix
+
+    for qubits in targets:
+        mermin_polynomial = get_mermin_polynomial(len(qubits))
+        readout_basis = get_readout_basis(mermin_polynomial)
+
+        for theta in thetas:
+            mermin_sequences = create_mermin_sequences(
+                platform, qubits, readout_basis=readout_basis, theta=theta
+            )
+            options = ExecutionParameters(nshots=params.nshots)
+            # TODO: use unrolling
+            for basis, sequence in mermin_sequences.items():
+                results = platform.execute_pulse_sequence(sequence, options=options)
+                frequencies = calculate_frequencies(results, qubits)
+                for state, frequency in enumerate(frequencies.values()):
+                    data.register_qubit(
+                        MerminType,
+                        tuple(qubits),
+                        dict(
+                            theta=np.array([theta]),
+                            basis=np.array([basis]),
+                            state=np.array([str(format(state, f"0{len(qubits)}b"))]),
+                            frequency=np.array([frequency]),
+                        ),
+                    )
+    return data
+
+
+def _fit(data: MerminData) -> MerminResults:
+    """Fitting for Mermin protocol."""
+    targets = data.targets
+    results = {qubits: [] for qubits in targets}
+    mitigated_results = {qubits: [] for qubits in targets}
+    basis = np.unique(data.data[targets[0]].basis)
+    for qubits in targets:
+        mermin_polynomial = get_mermin_polynomial(len(qubits))
+        mermin_coefficients = get_mermin_coefficients(mermin_polynomial)
+        for theta in data.thetas:
+            qubit_data = data.data[qubits]
+            outputs = []
+            mitigated_outputs = []
+            for base in basis:
+
+                data_filter = (qubit_data.basis == base) & (qubit_data.theta == theta)
+                state_freq = qubit_data[data_filter].frequency
+
+                outputs.append(
+                    {
+                        format(i, f"0{len(qubits)}b"): freq
+                        for i, freq in enumerate(state_freq)
+                    }
+                )
+
+                if data.mitigation_matrix:
+                    mitigated_output = np.dot(
+                        data.mitigation_matrix[qubits],
+                        state_freq,
+                    )
+                    mitigated_outputs.append(
+                        {
+                            format(i, f"0{len(qubits)}b"): freq
+                            for i, freq in enumerate(mitigated_output)
+                        }
+                    )
+            if data.mitigation_matrix:
+                mitigated_results[tuple(qubits)].append(
+                    compute_mermin(mitigated_outputs, mermin_coefficients)
+                )
+            results[tuple(qubits)].append(compute_mermin(outputs, mermin_coefficients))
+    return MerminResults(
+        mermin=results,
+        mermin_mitigated=mitigated_results,
+    )
+
+
+def _plot(data: MerminData, fit: MerminResults, target):
+    """Plotting function for Mermin protocol."""
+    figures = []
+    targets = data.targets
+
+    n_targets = len(targets)
+    classical_bound = 2 ** (n_targets // 2)
+    quantum_bound = 2 ** ((n_targets - 1) / 2) * (2 ** (n_targets // 2))
+
+    fig = go.Figure(layout_yaxis_range=[-quantum_bound - 0.2, quantum_bound + 0.2])
+    if fit is not None:
+        fig.add_trace(
+            go.Scatter(
+                x=data.thetas,
+                y=fit.mermin[tuple(target)],
+                name="Bare",
+            )
+        )
+        if fit.mermin_mitigated:
+            fig.add_trace(
+                go.Scatter(
+                    x=data.thetas,
+                    y=fit.mermin_mitigated[tuple(target)],
+                    name="Mitigated",
+                )
+            )
+
+    fig.add_trace(
+        go.Scatter(
+            mode="lines",
+            x=data.thetas,
+            y=[+classical_bound] * len(data.thetas),
+            line_color="gray",
+            name="Classical limit",
+            line_dash="dash",
+            legendgroup="classic",
+        )
+    )
+
+    fig.add_trace(
+        go.Scatter(
+            mode="lines",
+            x=data.thetas,
+            y=[-classical_bound] * len(data.thetas),
+            line_color="gray",
+            name="Classical limit",
+            legendgroup="classic",
+            line_dash="dash",
+            showlegend=False,
+        )
+    )
+
+    fig.add_trace(
+        go.Scatter(
+            mode="lines",
+            x=data.thetas,
+            y=[+quantum_bound] * len(data.thetas),
+            line_color="gray",
+            name="Quantum limit",
+            legendgroup="quantum",
+        )
+    )
+
+    fig.add_trace(
+        go.Scatter(
+            mode="lines",
+            x=data.thetas,
+            y=[-quantum_bound] * len(data.thetas),
+            line_color="gray",
+            name="Quantum limit",
+            legendgroup="quantum",
+            showlegend=False,
+        )
+    )
+
+    fig.update_layout(
+        xaxis_title="Theta [rad]",
+        yaxis_title="Mermin polynomial value",
+        xaxis=dict(range=[min(data.thetas), max(data.thetas)]),
+    )
+    figures.append(fig)
+
+    return figures, ""
+
+
+mermin = Routine(_acquisition, _fit, _plot)

src/qibocal/protocols/two_qubit_interaction/mermin/pulses.py

@@ -0,0 +1,84 @@
+from collections import defaultdict
+
+import numpy as np
+from qibolab.pulses import PulseSequence
+
+
+def create_mermin_sequence(platform, qubits, theta=None):
+    """Creates the pulse sequence to generate the bell states and with a theta-measurement"""
+
+    nqubits = len(qubits)
+    if theta is None:
+        theta = ((nqubits - 1) * 0.25 * np.pi) % (2 * np.pi)
+
+    virtual_z_phases = defaultdict(int)
+    sequence = PulseSequence()
+
+    for qubit in qubits:
+        sequence.add(
+            platform.create_RX90_pulse(
+                qubit, start=0, relative_phase=virtual_z_phases[qubit] + np.pi / 2
+            )
+        )
+
+    # TODO: Not hardcode topology
+
+    # qubits[0] needs to be the center qubit where everything is connected
+    for i in range(1, len(qubits)):
+        (cz_sequence1, cz_virtual_z_phases) = platform.create_CZ_pulse_sequence(
+            [qubits[0]] + [qubits[i]], sequence.finish + 8  # TODO: ask for the 8
+        )
+        sequence.add(cz_sequence1)
+        for qubit in cz_virtual_z_phases:
+            virtual_z_phases[qubit] += cz_virtual_z_phases[qubit]
+
+    t = sequence.finish + 8
+
+    for i in range(1, len(qubits)):
+        sequence.add(
+            platform.create_RX90_pulse(
+                qubits[i],
+                start=t,
+                relative_phase=virtual_z_phases[qubits[i]] - np.pi / 2,
+            )
+        )
+
+    virtual_z_phases[qubits[0]] -= theta
+
+    return sequence, virtual_z_phases
+
+
+def create_mermin_sequences(platform, qubits, readout_basis, theta):
+    """Creates the pulse sequences needed for the 4 measurement settings for chsh."""
+
+    mermin_sequences = {}
+
+    for basis in readout_basis:
+        sequence, virtual_z_phases = create_mermin_sequence(
+            platform, qubits, theta=theta
+        )
+        # t = sequence.finish
+        for i, base in enumerate(basis):
+            if base == "X":
+                sequence.add(
+                    platform.create_RX90_pulse(
+                        qubits[i],
+                        start=sequence.finish,
+                        relative_phase=virtual_z_phases[qubits[i]] + np.pi / 2,
+                    )
+                )
+            if base == "Y":
+                sequence.add(
+                    platform.create_RX90_pulse(
+                        qubits[i],
+                        start=sequence.finish,
+                        relative_phase=virtual_z_phases[qubits[i]],
+                    )
+                )
+        measurement_start = sequence.finish
+
+        for qubit in qubits:
+            sequence.add(platform.create_MZ_pulse(qubit, start=measurement_start))
+
+        mermin_sequences[basis] = sequence
+    return mermin_sequences