#!/usr/bin/env python
# coding: utf-8

# # Teleportation of a qubit across a quantum computer chip  
# This notebook teleports a qubit from one side of an IBM Eagle quantum computer chip to the other.  
# Some effort went into ensuring the code actually does this, and doesn't just teleport from one arbitrary qubit to another.

# In[1]:


import qiskit
import qiskit_ibm_runtime


# Get a quantum computer.  So far this is a simulator, to get the program working without eating real quantum computing resources.

# In[2]:


simulate = False
noise = True

if (noise or not simulate):
    import qiskit_ibm_runtime as ibm
    service = ibm.QiskitRuntimeService()
    ibm_backend = service.backend("ibm_kyoto") #pick this fairly noisy one to start.

if simulate:
    import qiskit_aer
    if not noise:
        backend = qiskit_aer.AerSimulator()
    else:
        backend = qiskit_aer.AerSimulator.from_backend(ibm_backend)
else:
    backend = ibm_backend


# In[3]:


num_qubits = 17 #17 is the most ibl_kyoto handles with reasonable noise.


# The coupling_map tells you the way the qubits are coupled for ecr gates.  e.g. [1,0] means you can do ecr with qubit 1 as the control and qubit 0 as the target.  If you do ecr in the other direction, it adds local operations and hence noise.

# In[4]:


if noise:
    coupling_map = ibm_backend.coupling_map
else:
    coupling_list = [(i, i+1) for i in range(num_qubits-1)]
    coupling_map = qiskit.transpiler.CouplingMap(coupling_list)

coupling_map.draw()


# Specify which physical qubits we will use for our logical qubits.  Choose them based on their noise, as e.g. Readout Assignment Error for a single qubit can be as high as 40%!

# In[5]:


#This is one path across the eagle chip, which seems not too noisy on a qubit and ecr level.  Check the latest calibration results on the chip you are using!
longest_layout = [13,12,17,30,29,28,35,47,46,45,54,64,63,62,72,81,80,79,91,98,97,96,109,114,113]

if noise: #We are based on the eagle chip, so can set the layout
    initial_layout = longest_layout[0:num_qubits]
else:
    initial_layout = range(num_qubits)

if len(initial_layout) != num_qubits:
    raise ValueError('Wrong number of qubits in initial layout')
#alternative method for specifying layout
#a, b = qubits
#initial_layout = Layout({a: 5, b: 6})


# Build a quantum circuit with n qubits and n classical bits (n odd).  The qubits are initialized as |0>'s.  We shall teleport from the first qubit to the last.

# In[6]:


qubits = qiskit.QuantumRegister(num_qubits, name = "qubits")
#One pair of cbits for each Bell Measurement
bell_cbits = [qiskit.ClassicalRegister(2, name = "bell_cbits_" + str(i)) for i in range(0, num_qubits-1, 2)]
final_cbit = qiskit.ClassicalRegister(1, name = "final_cbit")
qc = qiskit.QuantumCircuit(qubits, *bell_cbits, final_cbit)


# Setup a state to be teleported.  This is parameterized, so we can teleport many different states.  Thus far it's essentially $\frac{1}{\sqrt{2}}(|0> - i e^{i \phi} |1>)$.  Use native processor gates for the ibm eagle processor, to minimize the number of gates and hence the error.

# In[7]:


phi = qiskit.circuit.Parameter('phi')
qc.sx(0)
qc.rz(phi,0)


# Create a list of parameters for the initial state.

# In[8]:


import numpy
num_phis = 4
phi_values = numpy.linspace(0, numpy.pi/2, num_phis)


# Setup the pre-shared entanglement, using the native ecr direction.  We get the same entangled state whatever the native direction.

# In[9]:


def ecr_order(physical_qubits, edges):
    if physical_qubits in edges:
        entanglement_order = 0
    else:
        reversed_qubits = (physical_qubits[1], physical_qubits[0])
        if reversed_qubits in edges:
            entanglement_order = 1
        else:
            raise ValueError('These physical qubits are not linked')

    return entanglement_order


# In[10]:


edges = coupling_map.get_edges()

#Create the entangled pairs for logical qubits (1,2), (3,4) etc.
for i in range(1, qc.num_qubits, 2):
    physical_qubits = (initial_layout[i], initial_layout[i+1])
    entanglement_order = ecr_order(physical_qubits, edges)

    if entanglement_order == 0:
        qc.sx(i)
        qc.ecr(i,i+1)    
    else:
        qc.sx(i+1)
        qc.ecr(i+1,i)

qc.barrier()


# Display the circuit thus far: a state to be teleported and pre-shared entanglement.

# In[11]:


qc.draw(output='mpl')


# Perform one Bell Measurement on the initial qubit and the adjacent entangled qubit, and further ones to link the entangled pairs together, moving the state from qubit 0 to the final qubit.

# In[12]:


#Arrange the measurements to match the direction of the underlying hardware ecr gates.
measure_orders = []
for i in range(1, qc.num_qubits, 2):
    physical_qubits = (initial_layout[i-1], initial_layout[i])
    measure_order = ecr_order(physical_qubits, edges)
    measure_orders.append(measure_order) #We'll need these when interpreting the outcomes
    if measure_order == 0:
        qc.ecr(i-1,i)
        qc.sx(i-1)
    else:
        qc.ecr(i,i-1)
        qc.sx(i)
    
    qc.measure([i-1, i],[i-1,i])
qc.barrier()


# Draw the circuit so far:

# In[13]:


qc.draw(output='mpl')


# Next adjust the target qubit by the outcome of the Bell Measurement.  In the regular order, if the outcome is 00 or 10, do Z.  If outcome is 10 or 11, do X.  These operations commute, so if we do entanglement swapping we can easily add the bell measurement outcomes and apply X or I and Z or I to the final state: no need to do multiply X's, Z's or anything else. Similarly in the reverse order.

# In[14]:


target_num = num_qubits-1
num_pairs = round((num_qubits-1)/2)

import qiskit.circuit.classical.expr as expr

#Apply an X gate if x_expr is true, instead of many sequential X gates, potentially one for each pair.
#Similar for z gate.
for i in range(0, num_pairs):    
    c_register = bell_cbits[i]
    
    if measure_orders[i] == 0:
        pair_x_expr = c_register[0]
    else:        
        pair_x_expr = expr.equal(c_register[0], c_register[1]) #Flip if 00 or 11    

    pair_z_expr = expr.logic_not(c_register[1])
    
    if i==0:
        x_expr = pair_x_expr
        z_expr = pair_z_expr
    else:
        x_expr = expr.not_equal(pair_x_expr, x_expr) #two X gates make Identity.
        z_expr = expr.not_equal(pair_z_expr, z_expr)

with qc.if_test(x_expr) as else_:
    qc.x(target_num)
with qc.if_test(z_expr) as else_:
    qc.rz(numpy.pi, target_num)
qc.barrier()


# Finally measure $\sigma_y$ on the target qubit, i.e. in the |0>+i|1> vs |0>-i|1> basis.  This is chosen so that it changes with the input state, and can be done with a native .sx operation followed by measuring the computational basis.

# In[15]:


qc.sx(target_num)
qc.measure(target_num, target_num)


# Display our final logical circuit.

# In[16]:


qc.draw(output='mpl')


# Convert our circuit into native operations for the backend, specifying which logical qubit is assigned to which physical qubit.  To reduce noise, this should have as few layers as possible, and should look the same as our non-converted circuit as we attempted to use only native operations.

# In[17]:


qc_transpiled = qiskit.transpile(qc, backend, initial_layout=initial_layout)


# In[18]:


qc_transpiled.draw(output='mpl', idle_wires=False)


# Finally, we can run our experiment!

# In[19]:


load_results = True

if not load_results:
    #How do we set max_execution_time (defined in seconds)?  
    sampler = qiskit_ibm_runtime.SamplerV2(backend, options={"default_shots": 4096})
    job = sampler.run([(qc_transpiled, phi_values)])
    
    if not simulate:
        print(job.usage_estimation)


# The job might be in the queue for a few hours before it runs, so get the job_id so we can retreive it later.

# In[20]:


if not load_results:
    print(job.job_id())


# In case we need to load the job later

# In[20]:


if load_results:
    service = ibm.QiskitRuntimeService()
    job = service.job('crk9fvq14ys000888f00'); #The result of job_id() 


# In[21]:


results = job.result()


# Process the results.

# In[22]:


#get the measurement output counts just for the final (output) qubit
marginal_counts = [ results[0].data.final_cbit.get_counts(c) for c in range(results[0].data.final_cbit.size)]


# In[23]:


#One doc claimed you can't modify this when dynamic=True, but it does change as this shows.
num_shots = sum(marginal_counts[0].values())
print(num_shots)


# Calculate the expectation of the $\sigma_y$ operator.  $-1 \le E(\sigma_y) \le +1$.

# In[24]:


#To do sigma_y we applied sqrt(X) then measured $\sigma_z$.  Measurement outcome 0 corresponds to $\sigma_z = +1$, i.e. $\sigma_y = +1$.
E_sigma_y = [qiskit.result.sampled_expectation_value(counts, 'Z') for counts in marginal_counts]


# Plot the the output across various input states.

# In[25]:


import matplotlib
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(10, 6))

# results from hardware
ax.plot(phi_values / numpy.pi, E_sigma_y, "o-", label="QuantumSim", zorder=3)
#ax.plot(phis / numpy.pi, E_sigma_z, "o-", label="QuantumComp", zorder=3)

# classical bound +-0.5
ax.axhline(y=0.5, color="0.9", linestyle="--")
ax.axhline(y=-0.5, color="0.9", linestyle="--")

# quantum bound, +-1
ax.axhline(1, color="0.9", linestyle="-.")
ax.axhline(-1, color="0.9", linestyle="-.")
ax.fill_between(phi_values / numpy.pi, 1, 0.5, color="0.6", alpha=0.7)
ax.fill_between(phi_values / numpy.pi, -1, -0.5, color="0.6", alpha=0.7)

# set x tick labels to the unit of pi.
ax.xaxis.set_major_formatter(matplotlib.ticker.FormatStrFormatter('%g $\\pi$'))
ax.xaxis.set_major_locator(matplotlib.ticker.MultipleLocator(base=0.5))

# set labels, and legend
plt.xlabel("Phi")
plt.ylabel("Teleportation")
plt.legend()
plt.show()


# We are finished.  Note the quantum computer and library versions as they change often, and the code might not run on other versions.

# In[ ]:


backend.name


# In[ ]:


if not simulate:
    print(backend.processor_type)


# Note: the Eagle chip used in ibm_brisbane is approx 25mm wide.  
# We naively guess across the qubits inside the chip is perhaps 15mm.

# In[ ]:


qiskit_ibm_runtime.version.get_version_info()


# In[ ]:


qiskit.version.get_version_info()