Soluzione
# Imports
import numpy as np
import matplotlib.pyplot as plt
# Main qiskit imports
from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit_aer import AerSimulator
from qiskit_ibm_runtime.fake_provider import FakePerth
# Error mitigation
from qiskit_experiments.library.characterization import LocalReadoutError
# Noisy backend
backend = AerSimulator.from_backend(FakePerth())
# Local simulator
simulator = AerSimulator()
def zz_pump(q, c, p, system, ancilla):
"""Dà come risultato un QuantumCircuit che implementa il canale della pompa ZZ sui qubit di sistema
Args:
q (QuantumRegister): il registro da usare per il circuito
c (ClassicalRegister): il registro da usare per la misura dei qubit di sistema
p (float): l'efficienza del canale, fra 0 e 1
system (list): lista di indici per i qubit di sistema
ancilla (int): indice del qubit ausiliare
Returns:
Un oggetto QuantumCircuit
"""
zz = QuantumCircuit(q, c)
theta = 2 * np.arcsin(np.sqrt(p))
# Mappiamo l'informazione sugli ausiliari
zz.cx(q[system[0]], q[system[1]])
zz.x(q[ancilla])
zz.cx(q[system][1], q[ancilla])
# Rotazione condizionata
zz.cu(theta, 0.0, 0.0, 0.0, q[ancilla], q[system[1]])
# Mapping inverso
zz.cx(q[system[1]], q[ancilla])
# Misura
zz.h(q[system[0]])
zz.measure(q[system[0]], c[0])
zz.measure(q[system[1]], c[1])
return zz
def xx_pump(q, c, p, system, ancilla):
"""Dà come risultato un QuantumCircuit che implementa il canale della pompa ZZ sui qubit di sistema
Args:
q (QuantumRegister): il registro da usare per il circuito
c (ClassicalRegister): il registro da usare per la misura dei qubit di sistema
p (float): l'efficienza del canale, fra 0 e 1
system (list): lista di indici per i qubit di sistema
ancilla (int): indice del qubit ausiliare
Returns:
Un oggetto QuantumCircuit
"""
xx = QuantumCircuit(q, c)
theta = 2 * np.arcsin(np.sqrt(p))
# Mappiamo l'informazione sugli ausiliari
xx.cx(q[system[0]], q[system[1]])
xx.h(q[system[0]])
xx.x(q[ancilla])
xx.cx(q[system[0]], q[ancilla])
# Rotazione condizionata
xx.cu(theta, 0.0, 0.0, 0.0, q[ancilla], q[system[0]])
# Mapping inverso
xx.cx(q[system[0]], q[ancilla])
# Misura
xx.measure(q[system[0]], c[0])
xx.measure(q[system[1]], c[1])
return xx
def zz_xx_pump(q, c, p, system, ancillae):
"""Dà come risultato un QuantumCircuit che implementa il canale composto sui qubit di sistema
Args:
q (QuantumRegister): il registro da usare per il circuito
c (ClassicalRegister): il registro da usare per la misura dei qubit di sistema
p (float): l'efficienza del canale, fra 0 e 1
system (list): lista di indici per i qubit di sistema
ancilla (int): indice del qubit ausiliare
Returns:
Un oggetto QuantumCircuit
"""
zx = QuantumCircuit(q, c)
theta = 2 * np.arcsin(np.sqrt(p))
# pompa ZZ
# Mappiamo l'informazione sugli ausiliari
zx.cx(q[system[0]], q[system[1]])
zx.x(q[ancillae[0]])
zx.cx(q[system[1]], q[ancillae[0]])
# Rotazione condizionata
zx.cu(theta, 0.0, 0.0, 0.0, q[ancillae[0]], q[system[1]])
# Mapping inverso
zx.cx(q[system[1]], q[ancillae[0]])
# pompa XX
# Mappiamo l'informazione sugli ausiliari
zx.h(q[system[0]])
zx.x(q[ancillae[1]])
zx.cx(q[system[0]], q[ancillae[1]])
# Rotazione condizionata
zx.cu(theta, 0.0, 0.0, 0.0, q[ancillae[1]], q[system[0]])
# Mapping inverso
zx.cx(q[system[0]], q[ancillae[1]])
# Misura
zx.measure(q[system[0]], c[0])
zx.measure(q[system[1]], c[1])
return zx
Per semplicità definiamo una funzione che dà come risultato la preparazione dei quattro stati iniziali
def initial_conditions(q, system):
"""Dà come risultato un dizionario che contiene quattro oggetti QuantumCircuit che preparano il sistema a due qubit in stati iniziali diversi
Args:
q (QuantumRegister): il registro da usare per il circuito
system (list): lista di indici per i qubit di sistema
Returns:
Un dizionario con gli oggetti QuantumState dello stato iniziale e una lista di label
"""
# State labels
state_labels = ['00', '01', '10', '11']
ic = {}
for ic_label in state_labels:
ic[ic_label] = QuantumCircuit(q)
# |01>
ic['01'].x(q[system[0]])
# |10>
ic['10'].x(q[system[1]])
# |11>
ic['11'].x(q[system[0]])
ic['11'].x(q[system[1]])
return ic, state_labels
SHOTS = 8192
# I valori per p
p_values = np.linspace(0, 1, 10)
# Creiamo i circuiti quantistici
q = QuantumRegister(5, name='q')
c = ClassicalRegister(2, name='c')
## Indice del qubit di sistema
system = [2, 1]
## Indice dei qubit ausiliari
a_zz = 0
a_xx = 4
## Prepariamo i qubit in quattro condizioni iniziali
ic_circs, ic_state_labels = initial_conditions(q, system)
## Tre canali diversi, ciascuno con
## quattro condizioni iniziali e dieci valori di p
pumps = ['ZZ', 'XX', 'ZZ_XX']
circuits = {}
for pump in pumps:
circuits[pump] = {}
for ic in ic_state_labels:
circuits[pump][ic] = []
for ic in ic_state_labels:
for p in p_values:
circuits['ZZ'][ic].append(ic_circs[ic].compose(zz_pump(q, c, p, system, a_zz)))
circuits['XX'][ic].append(ic_circs[ic].compose(xx_pump(q, c, p, system, a_xx)))
circuits['ZZ_XX'][ic].append(ic_circs[ic].compose(zz_xx_pump(q, c, p, system, [a_zz, a_xx])))
circuits['ZZ_XX']['00'][1].draw(output='mpl')
# Eseguiamo i circuiti sul simulatore locale
jobs_sim = {}
for pump in pumps:
jobs_sim[pump] = {}
for ic in ic_state_labels:
jobs_sim[pump][ic] = simulator.run(circuits[pump][ic], shots = SHOTS)
# Analizziamo i risultati
overlaps_sim = {}
for pump in pumps:
overlaps_sim[pump] = {}
for ic in ic_state_labels:
overlaps_sim[pump][ic] = [0.0]*len(p_values)
for i in range(len(p_values)):
for ic in ic_state_labels:
counts = jobs_sim[pump][ic].result().get_counts(i)
for outcome in counts:
overlaps_sim[pump][outcome][i] += counts[outcome]/(4.0 * float(SHOTS))
# Rappresentiamo graficamente i risultati
fig_idx = 131
plt.figure(figsize=(15,6))
bell_labels = {'00': r"$| \phi^{+} \rangle$", '01': r"$| \phi^{-} \rangle$", '10': r"$| \psi^{+} \rangle$", '11': r"$| \psi^{-} \rangle$"}
for pump in pumps:
plt.subplot(fig_idx)
for outcome in overlaps_sim[pump]:
plt.plot(p_values, overlaps_sim[pump][outcome], label = bell_labels[outcome])
plt.xlabel('p')
plt.ylabel('Overlap')
fig_idx += 1
plt.grid()
plt.legend()
exp = LocalReadoutError(system, backend=backend)
exp.analysis.set_options(plot=True)
res = exp.run(shots=SHOTS)
mitigator = res.analysis_results(0).value
res.figure(0)
# Eseguiamo i circuiti sul simulatore locale
jobs_sim = {}
for pump in pumps:
jobs_sim[pump] = {}
for ic in ic_state_labels:
jobs_sim[pump][ic] = backend.run(circuits[pump][ic], shots = SHOTS)
overlaps = {}
overlaps_mit = {}
for pump in pumps:
overlaps[pump] = {}
overlaps_mit[pump] = {}
for ic in ic_state_labels:
overlaps[pump][ic] = [0.0]*len(p_values)
overlaps_mit[pump][ic] = [0.0]*len(p_values)
for i in range(len(p_values)):
for ic in ic_state_labels:
counts = jobs_sim[pump][ic].result().get_counts(i)
unmitigated_probs = {label: count / SHOTS for label, count in counts.items()}
mitigated_quasi_probs = mitigator.quasi_probabilities(unmitigated_probs)
mitigated_probs = mitigated_quasi_probs.nearest_probability_distribution().binary_probabilities()
for outcome in unmitigated_probs:
overlaps[pump][outcome][i] += unmitigated_probs[outcome]/4
for outcome in mitigated_probs:
overlaps_mit[pump][outcome[::-1]][i] += mitigated_probs[outcome]/4
# Plottiamo i risultati
fig_idx = 131
plt.figure(figsize=(15,6))
colors = plt.rcParams['axes.prop_cycle'][0:4]
cycler_2 = plt.cycler(color=colors)
bell_labels = {'00': r"$| \phi^{+} \rangle$", '01': r"$| \phi^{-} \rangle$", '10': r"$| \psi^{+} \rangle$", '11': r"$| \psi^{-} \rangle$"}
for pump in pumps:
plt.subplot(fig_idx)
plt.gca().set_prop_cycle(cycler_2)
for outcome in overlaps[pump]:
plt.plot(p_values, overlaps[pump][outcome], label = f"{bell_labels[outcome]} noisy", ls='--')
for outcome in overlaps_mit[pump]:
plt.plot(p_values, overlaps_mit[pump][outcome], label = f"{bell_labels[outcome]} mitigated")
plt.xlabel('p')
plt.ylabel('Overlap')
fig_idx += 1
plt.grid()
plt.legend()
La differenza tra i risultati con e senza mitigazione degli errori di lettura non è molto marcata, poiché in questo caso stiamo misurando solo 2 qubit. Più qubit abbiamo, più errori si accumulano e più diventa importante la mitigazione degli errori di lettura!