# If you already have saved IBM credentials with previous version, 
# update your credentials that is stored in disk
IBMQ.update_account()

# import
from qiskit import *
import numpy as np
from qiskit.visualization import plot_histogram, plot_state_city, plot_state_paulivec
from qiskit.quantum_info import Pauli, state_fidelity, basis_state, process_fidelity
from qiskit.tools.monitor import backend_overview, backend_monitor, job_monitor

# Create a quantum register with 3 qubits
q = QuantumRegister(3,'q')
# Form a quantum circuit
# Note that the circuit name is optional
qc = QuantumCircuit(q,name="first_qc")
# Display the quantum circuit
qc.draw()

# Create a classical register with 3 bits
c = ClassicalRegister(3,'c')

meas = QuantumCircuit(q,c,name="first_m")
meas.barrier(q)
meas.measure(q[0], c[0])
meas.measure(q[1], c[1])
meas.measure(q[2], c[2])
#meas.measure(q[3], c[3])
#meas.measure(q[4], c[4])
meas.draw()

# Add a x gate on qubit 0.
qc.x(0)
# Add a Hadamard gate on qubit 1, putting this in superposition.
qc.h(1)
# Add a CX (CNOT) gate on control qubit 1 
# and target qubit 2 to create an entangled state.
qc.cx(1, 2)
qc.draw()

# Quantum circuits can be added with + operations
# Add two pre-defined circuits
qc_comp=qc+meas
qc_comp.draw()

# Draw the quantum circuit in a different (slightly better) format
qc_comp.draw(output='mpl')

# Use Aer's qasm_simulator
backend_q = Aer.get_backend('qasm_simulator')

# Execute the circuit on the qasm simulator.
# We've set the number of repeats of the circuit
# to be 1024, which is default.
job_sim1 = execute(qc_comp, backend_q, shots=1024)

job_sim1.status()

# Grab the results from the job.
result_sim1 = job_sim1.result()
result_sim1

result_sim1.get_counts(qc_comp)

plot_histogram(result_sim1.get_counts(qc_comp))

# Use Aer's statevector_simulator
backend_sv = Aer.get_backend('statevector_simulator')

# Execute the circuit on the statevector simulator.
job_sim2 = execute(qc, backend_sv)

# Grab the results from the job.
result_sim2 = job_sim2.result()
# Output the entire result
job_sim2.result()

# See output state as a vector
outputstate = result_sim2.get_statevector(qc, decimals=5)
print(outputstate)

# Visualize density matrix
plot_state_city(outputstate)

# Examine state fidelity, compare with a basis state |111>
sf1=state_fidelity(basis_state('111', 3), outputstate)
# Examine state fidelity w.r.t the desired state
desired_state=[0,1,0,0,0,0,0,1]/np.sqrt(2)
sf2=state_fidelity(desired_state, outputstate)
print(sf1)
print(sf2)

# Create a simple quantum register with 2 qubit 
# to examine unitary simulator
qu2 = QuantumRegister(2,'qu')
# Form a quantum circuit
# Note that the circuit name is optional
qu = QuantumCircuit(qu2,name="unitary_qc")
# Add a single qubit rotation
qu.z(0)
qu.z(1)
# Display the quantum circuit
qu.draw(output='mpl')

# Use Aer's unitary_simulator
backend_u = Aer.get_backend('unitary_simulator')

# Execute the circuit on the unitary simulator.
job_usim = execute(qu, backend_u)

# Grab the results from the job.
result_usim = job_usim.result()
result_usim

# Output the unitary matrix
unitary = result_usim.get_unitary(qu)
print('%s\n' % unitary)
# Check process fidelity
pf1=process_fidelity(Pauli(label='IX').to_matrix(), unitary)
pf2=process_fidelity(Pauli(label='ZZ').to_matrix(), unitary)
print('%s\n' % pf1)
print('%s\n' % pf2)

# Measure an expectation value of a Pauli observable
Z = np.array([[1, 0], [0, -1]])
X = np.array([[0, 1], [1, 0]])
Y = np.array([[0, -1j], [1j, 0]])
ZZ=np.kron(Z,Z)
XX=np.kron(X,X)

# Or use qiskit's built-in function to define Pauli matrices
ZI = Pauli(label='ZI').to_matrix()
XI = Pauli(label='XI').to_matrix()
YI = Pauli(label='YI').to_matrix()

# Create some functions for measuring expectation values
# For pure states
def expectation_state(state, Operator):
    return np.dot(state.conj(), np.dot(Operator, state))
# For density matrices
def expectation_density(density, Operator):
    return np.trace(Operator @ density)

IBMQ.enable_account('YOUR_TOKERN')

print(IBMQ.stored_account())
print(IBMQ.active_account())

my_provider=IBMQ.get_provider()
my_provider.backends()

# Retrieve IBM Quantum device information
backend_overview()

# Retrieve a specific IBM Quantum device information
# ibmq_16_melbourne as an example
backend_mel = my_provider.get_backend('ibmq_16_melbourne')
backend_monitor(backend_mel)

# Create a 3-qubit GHZ state
ghz3 = QuantumRegister(3,'q')
qc= QuantumCircuit(ghz3)
qc.h(0)
qc.cx(0,1)
qc.cx(0,2)

qghz3=qc+meas
qghz3.draw(output='mpl')

# Define backend as one of the real devices (5-qubit systems for now)
backend_qx2 = my_provider.get_backend('ibmqx2')
backend_qx4 = my_provider.get_backend('ibmqx4')
# There is a new 5-qubit backend
backend_ou = my_provider.get_backend('ibmq_ourense')

# Run 3-qubit GHZ experiment on a 5-qubit device (try ourense)
job_exp = execute(qghz3, backend=backend_ou)
job_monitor(job_exp)

# Grab experimental results
result_ou = job_exp.result()
counts_ou = result_exp.get_counts(qghz3)
plot_histogram(counts_ou)

# Run the same experiment on the other 5-qubit device as well.
job_exp2 = execute(qghz3, backend=backend_qx4)
job_monitor(job_exp2)

# Grab experimental results
result_exp2 = job_exp2.result()
counts_exp2 = result_exp2.get_counts(qghz3)

# Finally, run the same experiment on the 14-qubit device.
job_exp3 = execute(qghz3, backend=backend_mel)
job_monitor(job_exp3)

# Grab experimental results
result_mel = job_exp3.result()
counts_mel = result_mel.get_counts(qghz3)
plot_histogram(counts_mel)

# Now, compare to theory by running it on qasm_simulator
job_qasm = execute(qghz3,backend=backend_q)
result_qasm = job_qasm.result()
counts_qasm = result_qasm.get_counts(qghz3)
# Plot both experimental and ideal results
plot_histogram([counts_qasm,counts_ou,counts_mel],
               color=['black','green','blue'],
               legend=['QASM','Ourense','Melbourne'])

# Create a quantum register with 3 qubits
q = QuantumRegister(3,'q')
# Form a quantum circuit
qent = QuantumCircuit(q)
qent.initialize([1, 0, 0, 0, 0, 0, 0, 1] / np.sqrt(2), [q[0], q[1], q[2]])
# Display the quantum circuit
qent.draw(output='mpl')

# Use Aer's state vector simulator
backend_sv = Aer.get_backend('statevector_simulator')

# Execute and plot output state
output_ent = execute(qent, backend_sv).result().get_statevector(0, decimals=3)
plot_state_city(output_ent)

# Build a sub-circuit
# sub_ghz(n,c,t) creates a quantum circuit for creating an n-qubit GHZ state.
# n is the total number of qubits, c is an index for a control qubit, t is a vector of indices for target qubits.
def sub_ghz(n,c,t):
    sub_q = QuantumRegister(n)
    sub_circ = QuantumCircuit(sub_q, name='sub_ghz')
    sub_circ.h(c)
    for i in t:
        sub_circ.cx(c,i)
    return sub_circ

# Create a sub-circuit for 3-qubit GHZ state prep.
sub_ghz3 = sub_ghz(3,0,[1,2])
sub_inst3 = sub_ghz3.to_instruction()

# Create a sub-circuit consisting of T gate and cnot.
sub_tcnot = QuantumCircuit(QuantumRegister(2),name='t+cnot')
sub_tcnot.t(0)
sub_tcnot.cx(0,1)

q = QuantumRegister(5, 'q')
circ = QuantumCircuit(q)

# Now create a quantum circuit using sub-circuits.
# Apply Hadamard on all gates.
for i in range(5):
    circ.h(i)
#
circ.barrier()
# Appy 2-qubit GHZ preparation on qubits 0 and 1, and on qubits 3 and 4.
circ.append(sub_tcnot, [q[0],q[1]])
circ.append(sub_tcnot, [q[3], q[4]])
#    
circ.barrier()
# Apply 3-qubit GHZ preparation on all 3 neighbouring qubits.
for i in range(3):
    circ.append(sub_inst3, [q[i], q[i+1],q[i+2]])
circ.draw(output='mpl')

# We can also decompose the circuit that is constructed using sub-circuits.
circ.decompose().draw(output='mpl')

from qiskit.circuit import Parameter

# Define two parameters, t1 and t2
theta1 = Parameter('t1')
theta2 = Parameter('t2')

# Build a 5-qubit circuit
qc = QuantumCircuit(5, 1)

# First parameter, t1, is used for a single qubit rotation of a controlled qubit
qc.ry(theta1,0)
for i in range(5-1):
    qc.cx(i, i+1)
qc.barrier()
# Second parameter, t2, is used here
qc.rz(theta2,range(5))
qc.barrier()

for i in reversed(range(5-1)):
    qc.cx(i, i+1)
    
qc.ry(theta1,0)
qc.measure(0, 0)

qc.draw(output='mpl')

theta1_range = np.linspace(0, np.pi/2, 2)
theta2_range = np.linspace(0, 2 * np.pi, 64)

circuits = [qc.bind_parameters({theta1: theta_val1, theta2: theta_val2})
            for theta_val1 in theta1_range for theta_val2 in theta2_range ]

# Visualize several circuits to check that correct circuits are generated correctly.
circuits[0].draw(output='mpl')
# circuits[1].draw(output='mpl')
# circuits[64].draw(output='mpl')

# theta1_range = np.linspace(0, np.pi/2, 2)
# theta2_range = np.linspace(0, 2 * np.pi, 64)

# circuits = [qc.bind_parameters({theta1: theta_val1, theta2: theta_val2})
#            for theta_val1 in theta1_range for theta_val2 in theta2_range ]

# Visualize several circuits to check that correct circuits are generated correctly.
# circuits[0].draw(output='mpl')
circuits[1].draw(output='mpl')
# circuits[64].draw(output='mpl')

# theta1_range = np.linspace(0, np.pi/2, 2)
# theta2_range = np.linspace(0, 2 * np.pi, 64)

# circuits = [qc.bind_parameters({theta1: theta_val1, theta2: theta_val2})
#            for theta_val1 in theta1_range for theta_val2 in theta2_range ]

# Visualize several circuits to check that correct circuits are generated correctly.
# circuits[0].draw(output='mpl')
# circuits[1].draw(output='mpl')
circuits[64].draw(output='mpl')

# Execute multiple circuits
job = execute(qc, backend=Aer.get_backend('qasm_simulator'),
              parameter_binds=[{theta1: theta_val1, theta2: theta_val2} 
                               for theta_val1 in theta1_range for theta_val2 in theta2_range])

# Store all counts
counts = [job.result().get_counts(i) for i in range(len(job.result().results))]

# Plot to visualize the result
import matplotlib.pyplot as plt
plt.figure(figsize=(12,6))
plt.plot(range(len(theta1_range)*len(theta2_range)),
         list(map(lambda counts: (counts.get('0',0)-counts.get('1',1))/1024,counts)))
plt.show()

# Reset operation forces the target qubit to go to |0>
q = QuantumRegister(1)
c = ClassicalRegister(1)
qc = QuantumCircuit(q, c)
qc.h(q)
qc.reset(q[0])
# qc.x(q)
qc.measure(q, c)
qc.draw(output='mpl')

# Test the result
job = execute(qc, backend=Aer.get_backend('qasm_simulator'), shots=1024)
job.result().get_counts(qc)

# Controlled qubit operation, controlled by classical bits
# This could be useful for designing a post-selection scheme of a quantum algorithm
q2 = QuantumRegister(2,'q')
cont = ClassicalRegister(1,'cont')
qc = QuantumCircuit(q2,cont)
qc.x(0)
qc.swap(q2[0],q2[1]).c_if(cont, 1)
qc.measure(q2[0],cont)
qc.draw(output='mpl')

job = execute(qc, backend=Aer.get_backend('qasm_simulator'), shots=1024)
job.result().get_counts(qc)

import random as r

qr = QuantumRegister(3)
cr1 = ClassicalRegister(1)
# Define a separate classical register for classically-controlled gate
cr2 = ClassicalRegister(1)
qtel = QuantumCircuit(qr,cr1,cr2)

# Make a reference circuit to compare
qref = QuantumCircuit(1,1)
# Prepare a 2-qubit bell state
qtel.h(1)
qtel.cx(1,2)
qtel.barrier()
# Prepare a random 1-qubit state
a1 = r.random()*np.pi
a2 = r.random()*np.pi
a3= r.random()*np.pi

qtel.u3(a1,a2,a3,0)

# Design a reference circuit for comparing to teleportation result
qref.u3(a1,a2,a3,0)
qref.measure(0,0)
print('Reference circuit:')
display(qref.draw(output='mpl'))

# Apply bell measurement on qubit 0 and 1
qtel.cx(0,1)
qtel.h(0)
qtel.barrier()
qtel.measure([0, 1],[0, 1])
qtel.z(2).c_if(cr1,1)
qtel.x(2).c_if(cr2,1)
qtel.measure(2,1)
print('Quantum teleportation circuit:')
display(qtel.draw(output='mpl'))

# n is the number of shots
n=10000
job_qtel = execute(qtel, backend_q, shots=n)
job_ref = execute(qref, backend_q, shots=n)

display(plot_histogram(job_ref.result().get_counts(qref)))
qtel_c = job_qtel.result().get_counts(qtel)
# Need some post-processing
print(qtel_c)
p1 = (qtel_c['1 0']+qtel_c['1 1'])/n
p0 = (qtel_c['0 0']+qtel_c['0 1'])/n

print('Prob. 1 from QT = %s' % p1)
print('Prob. 0 from QT = %s' % p0)

from qiskit.compiler import transpile
from qiskit.transpiler import PassManager, Layout
from qiskit.tools.visualization import plot_gate_map

# Let's take a look at the layout of physical devices with the gate directions
display(plot_gate_map(backend_qx4,plot_directed=True))
display(plot_gate_map(backend_ou,plot_directed=True))
display(plot_gate_map(backend_mel,plot_directed=True))

# Create a dummy circuit for default transpiler demonstration
qr = QuantumRegister(3)
cr = ClassicalRegister(3)
qc = QuantumCircuit(qr,cr)
qc.h(0)
qc.swap(0,1)
qc.cx(1,0)
qc.s(2)
qc.x(2)
qc.h(2)
qc.measure(qr[2],cr[2])
qc.h(0)
qc.cx(0,1)
qc.cx(0,1)
qc.h(0)
qc.measure(qr[0],cr[0])
qc.measure(qr[1],cr[1])
print('Original circuit')
display(qc.draw(output='mpl'))

# Transpile the circuit to run on ibmqx4
qt_qx4 = transpile(qc,backend_qx4)
print('Transpiled circuit for ibmqx4')
display(qt_qx4.draw(output='mpl'))

# Transpile the circuit to run on ibmq_Ourense
qt_ou = transpile(qc,backend_ou)
print('Transpiled circuit for ibmqx_Oursense')
display(qt_ou.draw(output='mpl'))

qr = QuantumRegister(3,'q')
qc = QuantumCircuit(qr)
qc.h(0)
qc.cx(0,1)
qc.cx(0,2)
qc.draw(output='mpl')

qc_t = transpile(qc, basis_gates=['u1','u2','u3','cx'])
# Display transpiled circuit
display(qc_t.draw(output='mpl'))

# Print out some circuit properties
# Total nunmber of operations
print("Number of operations = %s" % qc_t.size())
# Circuit depth
print("Circuit depth = %s" % qc_t.depth())
# Number of qubits
print("Number of qubits = %s" % qc_t.width())
# Breakdown of operations by type
print("Operations: %s" % qc_t.count_ops())
# Number of unentangled subcircuits in this circuit.
# In principle, each subcircuit can be executed on a different quantum device.
print("Number of unentangled subcircuits = %s" % qc_t.num_tensor_factors())

# Map it onto 5 qubit backend ibmqx2
qc_qx2 = transpile(qc, backend_qx2, basis_gates=['u1','u2','u3','cx'])
display(qc_qx2.draw(output='mpl'))
display(plot_gate_map(backend_qx2,plot_directed=True))

# Print out some circuit properties
# Total nunmber of operations
print("Number of operations = %s" % qc_qx2.size())
# Circuit depth
print("Circuit depth = %s" % qc_qx2.depth())

# Map it onto 5 qubit backend ibmqx4
qc_qx4 = transpile(qc, backend_qx4,basis_gates=['u1','u2','u3','cx'])
display(qc_qx4.draw(output='mpl'))
display(plot_gate_map(backend_qx4,plot_directed=True))

# Print out some circuit properties
# Total nunmber of operations
print("Number of operations = %s" % qc_qx4.size())
# Circuit depth
print("Circuit depth = %s" % qc_qx4.depth())

# Customize the layout
layout = Layout({qr[0]: 3, qr[1]: 4, qr[2]: 2})
# Map it onto 5 qubit backend ibmqx4
qc_qx4_new = transpile(qc, backend_qx4,initial_layout=layout,basis_gates=['u1','u2','u3','cx'])
display(qc_qx4_new.draw(output='mpl'))

# Print out some circuit properties
# Total nunmber of operations
print("Number of operations = %s" % qc_qx4_new.size())
# Circuit depth
print("Circuit depth = %s" % qc_qx4_new.depth())

# Test transpilation with Toffoli gate
qr = QuantumRegister(3,'q')
qccx = QuantumCircuit(qr)
qccx.ccx(0,1,2)
qccx.draw(output='mpl')

qccx_t = transpile(qccx, basis_gates=['u3','cx'])
# Display transpiled circuit
display(qccx_t.draw(output='mpl'))

# Print out some circuit properties
# Total nunmber of operations
print("Number of operations = %s" % qccx_t.size())
# Circuit depth
print("Circuit depth = %s" % qccx_t.depth())
# Number of qubits
print("Number of qubits = %s" % qccx_t.width())
# Breakdown of operations by type
print("Operations: %s" % qccx_t.count_ops())
# Number of unentangled subcircuits in this circuit.
# In principle, each subcircuit can be executed on a different quantum device.
print("Number of unentangled subcircuits = %s" % qccx_t.num_tensor_factors())

# Map it onto 5 qubit backend ibmqx2
qccx_qx2 = transpile(qccx, backend_qx2)
display(qccx_qx2.draw(output='mpl'))
display(plot_gate_map(backend_qx2,plot_directed=True))

# Print out some circuit properties
# Total nunmber of operations
print("Number of operations = %s" % qccx_qx2.size())
# Circuit depth
print("Circuit depth = %s" % qccx_qx2.depth())
# Number of cx gates
print("Number of cx operations = %s" % qccx_qx2.count_ops()['cx'])

# Map it onto 5 qubit backend ibmqx4
qccx_qx4 = transpile(qccx, backend_qx4)
display(qccx_qx4.draw(output='mpl'))
display(plot_gate_map(backend_qx4,plot_directed=True))

# Print out some circuit properties
# Total nunmber of operations
print("Number of operations = %s" % qccx_qx4.size())
# Circuit depth
print("Circuit depth = %s" % qccx_qx4.depth())
# Number of cx gates
print("Number of cx operations = %s" % qccx_qx4.count_ops()['cx'])

# Customize the layout
layout = Layout({qr[0]: 3, qr[1]: 4, qr[2]: 2})
# Map it onto 5 qubit backend ibmqx4
qccx_qx4_new = transpile(qccx, backend_qx4,initial_layout=layout)
display(qccx_qx4_new.draw(output='mpl'))
display(plot_gate_map(backend_qx4,plot_directed=True))

# Print out some circuit properties
# Total nunmber of operations
print("Number of operations = %s" % qccx_qx4_new.size())
# Circuit depth
print("Circuit depth = %s" % qccx_qx4_new.depth())
# Number of cx gates
print("Number of cx operations = %s" % qccx_qx4.count_ops()['cx'])

# The device information such as gate error and T1 and T2 also should be considered
# when customizing the layout.
backend_monitor(backend_qx2)

display(plot_gate_map(backend_ou,plot_directed=True))

# Customize the layout IBMQ Ourense
layout_ou = Layout({qr[0]: 2, qr[1]: 3, qr[2]: 1})
# Map it onto 5 qubit backend ibmqx4
qccx_ou = transpile(qccx, backend_ou,initial_layout=layout_ou)
display(qccx_ou.draw(output='mpl'))

# Print out some circuit properties
# Total nunmber of operations
print("Number of operations = %s" % qccx_ou.size())
# Circuit depth
print("Circuit depth = %s" % qccx_ou.depth())
# Number of cx gates
print("Number of cx operations = %s" % qccx_ou.count_ops()['cx'])

# Apply Toffoli among 5 qubits
qr = QuantumRegister(5,'q')
toff = QuantumCircuit(qr)

for i in range(5):
    toff.h(i)

toff.ccx(0,1,2)
toff.ccx(1,2,3)
toff.ccx(2,3,4)
    
display(toff.draw(output='mpl'))

# Print out some circuit properties
# Total nunmber of operations
print("Number of operations = %s" % toff.size())
# Count different types of operations
print("Operation counts = %s" % toff.count_ops())
# Circuit depth
print("Circuit depth = %s" % toff.depth())

toff_t = transpile(toff, backend_mel)
display(toff_t.draw(output='latex'))

# Print out some circuit properties
# Total nunmber of operations
print("Number of operations = %s" % toff_t.size())
# Count different types of operations
print("Operation counts = %s" % toff_t.count_ops())
# Circuit depth
print("Circuit depth = %s" % toff_t.depth())

# Transpile many times (100 times in this example) and pick the best one
tcircs0 = transpile([toff]*100, backend_mel,optimization_level=0)
tcircs1 = transpile([toff]*100, backend_mel,optimization_level=1)
tcircs2 = transpile([toff]*100, backend_mel,optimization_level=2)
tcircs3 = transpile([toff]*100, backend_mel,optimization_level=3)

num_cx0 = [c.count_ops()['cx'] for c in tcircs0]
num_cx1 = [c.count_ops()['cx'] for c in tcircs1]
num_cx2 = [c.count_ops()['cx'] for c in tcircs2]
num_cx3 = [c.count_ops()['cx'] for c in tcircs3]

num_tot0 = [c.size() for c in tcircs0]
num_tot1 = [c.size() for c in tcircs1]
num_tot2 = [c.size() for c in tcircs2]
num_tot3 = [c.size() for c in tcircs3]

num_depth0 = [c.depth() for c in tcircs0]
num_depth1 = [c.depth() for c in tcircs1]
num_depth2 = [c.depth() for c in tcircs2]
num_depth3 = [c.depth() for c in tcircs3]

plt.rcParams.update({'font.size': 16})

# Plot the number of CNOT gates
plt.figure(figsize=(12,6))
plt.plot(range(len(num_cx0)),num_cx0,'r',label='level 0')
plt.plot(range(len(num_cx1)),num_cx1,'b',label='level 1')
plt.plot(range(len(num_cx2)),num_cx2,'g',label='level 2')
plt.plot(range(len(num_cx3)),num_cx3,'k',label='level 3')
plt.legend(loc='upper left')
plt.xlabel('Random trial')
plt.ylabel('# of cx gates')
plt.show()

# Plot total number of gates
plt.figure(figsize=(12,6))
plt.plot(range(len(num_cx0)),num_tot0,'r',label='level 0')
plt.plot(range(len(num_cx1)),num_tot1,'b',label='level 1')
plt.plot(range(len(num_cx2)),num_tot2,'g',label='level 2')
plt.plot(range(len(num_cx3)),num_tot3,'k',label='level 3')
plt.legend(loc='upper left')
plt.xlabel('Random trial')
plt.ylabel('# of total gates')
plt.show()

# Plot the number of CNOT gates
plt.figure(figsize=(12,6))
plt.plot(range(len(num_cx0)),num_depth0,'r',label='level 0')
plt.plot(range(len(num_cx1)),num_depth1,'b',label='level 1')
plt.plot(range(len(num_cx2)),num_depth2,'g',label='level 2')
plt.plot(range(len(num_cx3)),num_depth3,'k',label='level 3')
plt.legend(loc='upper left')
plt.xlabel('Random trial')
plt.ylabel('Circuit depth')
plt.show()

print('Opt0: Minimum # of cx gates = %s' % min(num_cx0))
print('Opt0: The best circuit is the circut %s' % num_cx0.index(min(num_cx0)))

print('Opt1: Minimum # of cx gates = %s' % min(num_cx1))
print('Opt1: The best circuit is the circut %s' % num_cx1.index(min(num_cx1)))

print('Opt2: Minimum # of cx gates = %s' % min(num_cx2))
print('Opt2: The best circuit is the circut %s' % num_cx2.index(min(num_cx2)))

print('Opt3: Minimum # of cx gates = %s' % min(num_cx3))
print('Opt3: The best circuit is the circut %s' % num_cx3.index(min(num_cx3)))

# Create a 10-qubit GHZ state
qr = QuantumRegister(10,'q')
ghz = QuantumCircuit(qr)
ghz.h(0)
for i in range(9):
    ghz.cx(i,i+1)

display(ghz.draw(output='mpl'))

# Print out some circuit properties
# Total nunmber of operations
print("Number of operations = %s" % ghz.size())
# Count different types of operations
print("Operation counts = %s" % ghz.count_ops())
# Circuit depth
print("Circuit depth = %s" % ghz.depth())

ghz_t = transpile(ghz, backend_mel)
ghz_t.draw(output='latex')

# Print out some circuit properties
# Total nunmber of operations
print("Number of operations = %s" % ghz_t.size())
# Count different types of operations
print("Operation counts = %s" % ghz_t.count_ops())
# Circuit depth
print("Circuit depth = %s" % ghz_t.depth())

# Transpile 100 times and pick the best one
circs = transpile([ghz]*100, backend_mel)
num_cx = [c.count_ops()['cx'] for c in circs]
# Plot the number of CNOT gates
plt.figure(figsize=(12,6))
plt.plot(range(len(num_cx)),num_cx)
plt.show()
print('Minimum # of cx gates = %s' % min(num_cx))
print('The best circuit is the circut %s' % num_cx.index(min(num_cx)))

