Source code for qlauncher.routines.qiskit.mitigation_suppression.mitigation

  1from __future__ import annotations
  2
  3from collections import defaultdict
  4from itertools import chain
  5from typing import TYPE_CHECKING, Literal
  6
  7import numpy as np
  8from qiskit import QuantumCircuit, transpile
  9from qiskit._accelerate.circuit import CircuitInstruction as AccelerateInstruction
 10from qiskit.circuit import Instruction, Operation
 11
 12from qlauncher.routines.circuits import CIRCUIT_FORMATS
 13from qlauncher.utils import sum_counts
 14
 15from .base import CircuitExecutionMethod
 16
 17if TYPE_CHECKING:
 18    from qiskit.quantum_info import SparsePauliOp
 19
 20    from qlauncher.routines.qiskit.adapters import GateCircuitBackend
 21    from qlauncher.routines.qiskit.backends.qiskit_backend import QiskitBackend
 22
 23


[docs]
 24class NoMitigation(CircuitExecutionMethod):
 25    compatible_circuit = CIRCUIT_FORMATS
 26


[docs]
 27    def sample(self, circuit: CIRCUIT_FORMATS, backend: GateCircuitBackend, shots: int = 1024) -> dict[str, int]:
 28        return backend.sampler.run([circuit], shots=shots).result()[0].join_data().get_counts()


 29


[docs]
 30    def estimate(self, circuit: QuantumCircuit, observable: SparsePauliOp, backend: QiskitBackend) -> float:
 31        return backend.estimator.run([(circuit, observable)]).result()[0].data.evs




 32
 33


[docs]
 34class WeighedMitigation(CircuitExecutionMethod):
 35    def __init__(self, mitigation_methods: list[CircuitExecutionMethod], method_weights: list[float] | None = None) -> None:
 36        if method_weights is None:
 37            method_weights = [1.0] * len(mitigation_methods)
 38        if len(method_weights) != len(mitigation_methods):
 39            raise ValueError(
 40                f'You must provide as many weights as there are methods! Expected {len(mitigation_methods)}, got {len(method_weights)}'
 41            )
 42        self.weights = method_weights
 43        self.methods = mitigation_methods
 44        # TODO: Change into intersection between all mitigation methods
 45        self.compatible_circuit = mitigation_methods[0].compatible_circuit
 46        super().__init__()
 47


[docs]
 48    def sample(self, circuit: QuantumCircuit, backend: QiskitBackend, shots: int = 1024) -> dict[str, int]:
 49        counts = (m.sample(circuit, backend, shots) for m in self.methods)
 50        result = defaultdict(int)
 51        for count, weight in zip(counts, self.weights, strict=True):
 52            for k, v in count.items():
 53                result[k] += int(round(v * weight, 0))
 54        return result


 55


[docs]
 56    def estimate(self, circuit: QuantumCircuit, observable: SparsePauliOp, backend: QiskitBackend) -> float:
 57        return sum(m.estimate(circuit, observable, backend) * w for m, w in zip(self.methods, self.weights, strict=True)) / len(
 58            self.weights
 59        )




 60
 61


[docs]
 62class PauliTwirling(CircuitExecutionMethod):
 63    """
 64    Error mitigation technique based on averaging the results of running multiple "twirled" versions of the initial circuit.
 65    The method appends additional gates on both sides of random 2 qubit gates (cx, ecr).
 66    """
 67
 68    compatible_circuit = QuantumCircuit
 69
 70    def __init__(self, num_random_circuits: int, max_substitute_gates_per_circuit: int = 4, do_transpile: bool = True) -> None:
 71        self.num_random_circuits = num_random_circuits
 72        self.max_substitute_gates_per_circuit = max_substitute_gates_per_circuit
 73        self.do_transpile = do_transpile
 74
 75    def _random_replacement_op(self, inst: AccelerateInstruction) -> list[AccelerateInstruction]:
 76        op: Operation = inst.operation
 77        match op.name:
 78            case 'cx':
 79                return [
 80                    [
 81                        AccelerateInstruction(
 82                            operation=Instruction(name='x', num_qubits=1, num_clbits=0, params=[]), qubits=[inst.qubits[0]]
 83                        ),
 84                        AccelerateInstruction(
 85                            operation=Instruction(name='y', num_qubits=1, num_clbits=0, params=[]), qubits=[inst.qubits[1]]
 86                        ),
 87                        inst,
 88                        AccelerateInstruction(
 89                            operation=Instruction(name='y', num_qubits=1, num_clbits=0, params=[]), qubits=[inst.qubits[0]]
 90                        ),
 91                        AccelerateInstruction(
 92                            operation=Instruction(name='z', num_qubits=1, num_clbits=0, params=[]), qubits=[inst.qubits[1]]
 93                        ),
 94                    ],
 95                    [
 96                        AccelerateInstruction(
 97                            operation=Instruction(name='x', num_qubits=1, num_clbits=0, params=[]), qubits=[inst.qubits[0]]
 98                        ),
 99                        inst,
100                        AccelerateInstruction(
101                            operation=Instruction(name='x', num_qubits=1, num_clbits=0, params=[]), qubits=[inst.qubits[0]]
102                        ),
103                        AccelerateInstruction(
104                            operation=Instruction(name='x', num_qubits=1, num_clbits=0, params=[]), qubits=[inst.qubits[1]]
105                        ),
106                    ],
107                    [
108                        AccelerateInstruction(
109                            operation=Instruction(name='z', num_qubits=1, num_clbits=0, params=[]), qubits=[inst.qubits[0]]
110                        ),
111                        inst,
112                        AccelerateInstruction(
113                            operation=Instruction(name='z', num_qubits=1, num_clbits=0, params=[]), qubits=[inst.qubits[0]]
114                        ),
115                    ],
116                ][int(np.random.default_rng().integers(0, 3))]
117
118            case 'ecr':
119                return [
120                    [
121                        AccelerateInstruction(
122                            operation=Instruction(name='x', num_qubits=1, num_clbits=0, params=[]), qubits=[inst.qubits[0]]
123                        ),
124                        AccelerateInstruction(
125                            operation=Instruction(name='y', num_qubits=1, num_clbits=0, params=[]), qubits=[inst.qubits[1]]
126                        ),
127                        inst,
128                        AccelerateInstruction(
129                            operation=Instruction(name='x', num_qubits=1, num_clbits=0, params=[]), qubits=[inst.qubits[0]]
130                        ),
131                        AccelerateInstruction(
132                            operation=Instruction(name='y', num_qubits=1, num_clbits=0, params=[]), qubits=[inst.qubits[1]]
133                        ),
134                    ],
135                    [
136                        AccelerateInstruction(
137                            operation=Instruction(name='x', num_qubits=1, num_clbits=0, params=[]), qubits=[inst.qubits[0]]
138                        ),
139                        AccelerateInstruction(
140                            operation=Instruction(name='z', num_qubits=1, num_clbits=0, params=[]), qubits=[inst.qubits[1]]
141                        ),
142                        inst,
143                        AccelerateInstruction(
144                            operation=Instruction(name='x', num_qubits=1, num_clbits=0, params=[]), qubits=[inst.qubits[0]]
145                        ),
146                        AccelerateInstruction(
147                            operation=Instruction(name='z', num_qubits=1, num_clbits=0, params=[]), qubits=[inst.qubits[1]]
148                        ),
149                    ],
150                    [
151                        AccelerateInstruction(
152                            operation=Instruction(name='x', num_qubits=1, num_clbits=0, params=[]), qubits=[inst.qubits[1]]
153                        ),
154                        inst,
155                        AccelerateInstruction(
156                            operation=Instruction(name='x', num_qubits=1, num_clbits=0, params=[]), qubits=[inst.qubits[1]]
157                        ),
158                    ],
159                ][int(np.random.default_rng().integers(0, 3))]
160            case _:
161                return [inst]
162
163    def _twirl_circuit(self, transpiled_circuit: QuantumCircuit) -> QuantumCircuit:
164        """Apply random self.max_substitute_gates_per_circuit twirls on random (no replacement) gates of the circuit."""
165        circuit = transpiled_circuit.copy()
166
167        double_gates_with_indices: list[tuple[int, AccelerateInstruction]] = [
168            (i, x) for i, x in enumerate(circuit.data) if x.operation.num_qubits == 2
169        ]
170
171        choice_idxs = np.random.default_rng().choice(
172            range(len(double_gates_with_indices)),
173            size=min(self.max_substitute_gates_per_circuit, len(double_gates_with_indices)),
174            replace=False,
175        )
176
177        data_cpy = [[x] for x in circuit.data]
178
179        for i in choice_idxs:
180            data_cpy[i] = self._random_replacement_op(data_cpy[i][0])
181
182        circuit.data = list(chain.from_iterable(data_cpy))  # Collapse [[e1],[e2,e3],[e4],...] to [e1,e2,e3,e4,...]
183
184        return circuit
185
186    def _get_workable_circuit(self, circuit: QuantumCircuit, backend: QiskitBackend) -> QuantumCircuit:
187        """Get either transpiled circuit if do_transpile is set or copy of circuit you can change as you wish."""
188        return transpile(circuit, basis_gates=list(backend.backendv1v2.target.operation_names)) if self.do_transpile else circuit.copy()
189


[docs]
190    def sample(self, circuit: QuantumCircuit, backend: QiskitBackend, shots: int = 1024) -> dict[str, int]:
191        input_circ = self._get_workable_circuit(circuit, backend)
192        results = backend.sampler.run(
193            [transpile(self._twirl_circuit(input_circ), backend.backendv1v2) for _ in range(self.num_random_circuits)],
194            shots=shots // self.num_random_circuits,
195        ).result()
196
197        counts = [r.join_data().get_counts() for r in results]
198
199        return sum_counts(*counts)


200


[docs]
201    def estimate(self, circuit: QuantumCircuit, observable: SparsePauliOp, backend: QiskitBackend) -> float:
202        input_circ = self._get_workable_circuit(circuit, backend)
203
204        results = backend.estimator.run(
205            [
206                (transpile(self._twirl_circuit(input_circ), basis_gates=list(backend.backendv1v2.target.operation_names)), observable)
207                for _ in range(self.num_random_circuits)
208            ],
209        ).result()
210
211        sum_evs = 0
212        for r in results:
213            sum_evs += r.data.evs
214
215        return sum_evs / self.num_random_circuits




216
217


[docs]
218class ZeroNoiseExtrapolation(CircuitExecutionMethod):
219    """
220    Error mitigation technique based on fitting a model to data generated
221    by running a circuit made to multiply the error of the original circuit,
222    then predicting the values at x=0 (original circuit)
223    """
224
225    compatible_circuit = QuantumCircuit
226
227    def __init__(self, num_extrapolations: int = 4, polynomial_degree: int = 3, mode: Literal['linear', 'exponential'] = 'linear') -> None:
228        """
229        Args:
230                num_extrapolations (int, optional): Number of times the whole circuit is repeated for the largest X. Defaults to 4.
231                polynomial_degree (int, optional): Degree of fitted polynomial. Defaults to 3.
232                mode (Literal["linear", "exponential"], optional):
233                        Scaling method. "linear" keeps the original values as is,
234                        "exponential" applies log before fitting model then applies exp to the model prediction.
235                        Defaults to "linear".
236
237        Raises:
238                ValueError: If the polynomial degree is larger or equal to the number of data points (num_extrapolations)
239        """
240        super().__init__()
241        if polynomial_degree >= num_extrapolations:
242            raise ValueError('Degree must be lower than number of data points.')
243        self.num_extrapolations = num_extrapolations
244        self.degree = polynomial_degree
245        self.mode: Literal['linear', 'exponential'] = mode
246
247    def _get_repeated_circuits(self, circuit: QuantumCircuit) -> list[QuantumCircuit]:
248        result = []
249        meas_circ = circuit.copy()
250
251        mod_circuit: QuantumCircuit = circuit.remove_final_measurements(inplace=False)
252        inv = mod_circuit.inverse(annotated=True)
253
254        for _ in range(1, self.num_extrapolations + 1):
255            meas_circ.compose(inv, front=True, inplace=True)
256            meas_circ.compose(mod_circuit, front=True, inplace=True)
257            result.append(meas_circ.copy())
258
259        return result
260
261    def _get_zero_estimate(self, y_values: np.ndarray) -> float:
262        return np.polynomial.Polynomial.fit(np.array(range(1, self.num_extrapolations + 1)), y_values, self.degree).convert().coef[0]
263
264    def _get_zero_estimate_sampling(self, y_values: np.ndarray) -> np.ndarray:
265        return np.array([self._get_zero_estimate(y_values.T[i]) for i in range(len(y_values[0]))])
266
267    def _get_np_array_from_counts_dict(self, int_counts: dict[int, int], num_measured: int) -> np.ndarray:
268        result = np.zeros(2**num_measured)
269        for value, counts in int_counts.items():
270            result[value] = counts
271        return result
272


[docs]
273    def sample(self, circuit: QuantumCircuit, backend: GateCircuitBackend, shots: int = 1024) -> dict[str, int]:
274        counts = np.array(
275            [
276                self._get_np_array_from_counts_dict(res.join_data().get_int_counts(), circuit.num_clbits)
277                for res in backend.sampler.run(self._get_repeated_circuits(circuit), shots=shots).result()
278            ]
279        )
280
281        if self.mode == 'exponential':
282            counts = np.log(counts)
283            counts_fit = np.exp(self._get_zero_estimate_sampling(counts))
284        else:
285            counts_fit = np.maximum(self._get_zero_estimate_sampling(counts), 0) + (10 if counts.sum() == 0 else 0)
286
287        counts_dict = {}
288        for i, val in enumerate(counts_fit):
289            counts_dict[np.binary_repr(i, circuit.num_clbits)] = int(val)
290
291        return counts_dict


292


[docs]
293    def estimate(self, circuit: QuantumCircuit, observable: SparsePauliOp, backend: QiskitBackend) -> float:
294        evs = np.array(
295            [res.data.evs for res in backend.estimator.run([(x, observable) for x in self._get_repeated_circuits(circuit)]).result()]
296        )
297        if self.mode == 'exponential':
298            evs = np.log(evs)
299            return np.exp(self._get_zero_estimate(evs))
300        return self._get_zero_estimate(evs)