# useful additional packages 
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
import time
from pprint import pprint

# importing Qiskit
from qiskit import Aer, IBMQ
from qiskit.backends.ibmq import least_busy
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister, execute

# import basic plot tools
from qiskit.tools.visualization import plot_histogram

IBMQ.load_accounts()

"""
Explanation: <img src="../../../images/qiskit-heading.gif" alt="Note: In order for images to show up in this jupyter notebook you need to select File => Trusted Notebook" width="500 px" align="left">
W state in multi-qubit systems
The latest version of this notebook is available on https://github.com/qiskit/qiskit-tutorial.
For more information about how to use the IBM Q experience (QX), consult the tutorials, or check out the community.

Contributors
Pierre Decoodt, Université Libre de Bruxelles
End of explanation
"""


"Choice of the backend"
# using local qasm simulator
backend = Aer.get_backend('qasm_simulator')  

# using IBMQ qasm simulator 
# backend = IBMQ.get_backend('ibmq_qasm_simulator')
# using real device
# backend = least_busy(IBMQ.backends(simulator=False))

flag_qx2 = True
if backend.name() == 'ibmqx4':
        flag_qx2 = False
        
print("Your choice for the backend is: ", backend, "flag_qx2 is: ", flag_qx2)

# Here, two useful routine
# Define a F_gate
def F_gate(circ,q,i,j,n,k) :
    theta = np.arccos(np.sqrt(1/(n-k+1)))
    circ.ry(-theta,q[j])       
    circ.cz(q[i],q[j])
    circ.ry(theta,q[j])
    circ.barrier(q[i])
# Define the cxrv gate which uses reverse CNOT instead of CNOT
def  cxrv(circ,q,i,j) :
    circ.h(q[i])
    circ.h(q[j])
    circ.cx(q[j],q[i])
    circ.h(q[i])
    circ.h(q[j])
    circ.barrier(q[i],q[j])

"""
Explanation: Theoretical background
In addition to the GHZ states, the generalized W states, as proposed by Dür, Vidal and Cirac, in 2000, is a class of interesting examples of multiple qubit entanglement.
A generalized $n$ qubit W state can be written as :
$$ |W_{n}\rangle \; = \; \sqrt{\frac{1}{n}} \: (\:|10...0\rangle \: +  |01...0\rangle \: +...+ |00...1\rangle \:) $$
Here are presented circuits allowing to deterministically produce respectively a three, a four and a five qubit W state.
A 2016 paper by Firat Diker proposes an algorithm in the form of nested boxes allowing the deterministic construction of W states of any size $n$. The experimental setup proposed by the author is essentially an optical assembly including half-wave plates. The setup includes $n-1$ so-called two-qubit $F$ gates (not to be confounded with the Fredkin's three-qubit gate).
It is possible to construct the equivalent of such a $F$ gate on a superconducting quantum computing system using transmon qubits in ground and excited states. A $F_{k,\, k+1}$ gate with control qubit $q_{k}$ and target qubit $q_{k+1}$ is obtained here by:

First a rotation round Y-axis $R_{y}(-\theta_{k})$ applied on $q_{k+1}$
Then a controlled Z-gate $cZ$ in any direction between the two qubits $q_{k}$ and $q_{k+1}$
Finally a rotation round Y-axis $R_{y}(\theta_{k})$ applied on $q_{k+1}$

The matrix representations of a $R_{y}(\theta)$ rotation and of the $cZ$ gate can be found in the "Quantum gates and linear algebra" Jupyter notebook of the Qiskit tutorial. 
The value of $\theta_{k}$ depends on $n$ and $k$ following the relationship:
$$\theta_{k} = \arccos \left(\sqrt{\frac{1}{n-k+1}}\right) $$
Note that this formula for $\theta$ is different from the one mentioned in the Diker's paper. This is due to the fact that we use here Y-axis rotation matrices instead of $W$ optical gates composed of half-wave plates.
At the beginning, the qubits are placed in the state:  $|\varphi_{0} \rangle \, = \, |10...0 \rangle$.
This is followed by the application of $n-1$ sucessive $F$ gates. 
$$|\varphi_{1}\rangle = F_{n-1,\,n}\, ... \, F_{k,\, k+1}\, ... \, F_{2,\, 3}  \,F_{1,\, 2}\,|\varphi_{0} \rangle \,= \; \sqrt{\frac{1}{n}} \: (\:|10...0\rangle \: +  |11...0\rangle \: +...+ |11...1\rangle \:) $$
Then, $n-1$ $cNOT$ gates are applied. The final circuit is: 
$$|W_{n}\rangle \,= cNOT_{n,\, n-1}\, cNOT_{n-1,\, n-2}...cNOT_{k,\, k-1}...cNOT_{2,\, 1}\,\,|\varphi_{1} \rangle$$
Let's launch now in the adventure of producing deterministically W states, on simulator or in the real world!
Now you will have the opportunity to choose your backend.
(If you run the following cells in sequence, you will end with the local simulator, which is a good choice for a first trial).
End of explanation
"""


# 3-qubit W state Step 1
n = 3
q = QuantumRegister(n)
c = ClassicalRegister(n)
W_states = QuantumCircuit(q,c) 

W_states.x(q[2]) #start is |100>
F_gate(W_states,q,2,1,3,1) # Applying F12
    
for i in range(3) :
    W_states.measure(q[i] , c[i]) 

# circuits = ['W_states']

shots = 1024
time_exp = time.strftime('%d/%m/%Y %H:%M:%S')
print('start W state 3-qubit (step 1) on', backend, "N=", shots,time_exp)
result = execute(W_states, backend=backend, shots=shots)
time_exp = time.strftime('%d/%m/%Y %H:%M:%S')
print('end   W state 3-qubit (step 1) on', backend, "N=", shots,time_exp)
plot_histogram(result.result().get_counts(W_states))

"""
Explanation: Three-qubit W state, step 1
In this section, the production of a three qubit W state will be examined step by step.
In this circuit, the starting state is now: $ |\varphi_{0} \rangle \, = \, |100\rangle$.
The entire circuit corresponds to:                                 
$$ |W_{3}\rangle \,=\,  cNOT_{3,2}\,  \, cNOT_{2,1}\,  \, F_{2,3} \, 
\, F_{1,2} \,  \, |\varphi_{0} \rangle \, $$    
Run the following cell to see what happens when we first apply $F_{1,2}$.
End of explanation
"""


# 3-qubit W state, first and second steps
n = 3
q = QuantumRegister(n) 
c = ClassicalRegister(n)
W_states = QuantumCircuit(q,c) 
        
W_states.x(q[2]) #start is |100>
F_gate(W_states,q,2,1,3,1) # Applying F12
F_gate(W_states,q,1,0,3,2) # Applying F23
       
for i in range(3) :
    W_states.measure(q[i] , c[i]) 

shots = 1024
time_exp = time.strftime('%d/%m/%Y %H:%M:%S')
print('start W state 3-qubit (steps 1 + 2) on', backend, "N=", shots,time_exp)
result = execute(W_states, backend=backend, shots=shots)
time_exp = time.strftime('%d/%m/%Y %H:%M:%S')
print('end   W state 3-qubit (steps 1 + 2) on', backend, "N=", shots,time_exp)
plot_histogram(result.result().get_counts(W_states))

"""
Explanation: Three-qubit W state: adding step 2
In the previous step you obtained an histogram compatible with the following state:
$$ |\varphi_{1} \rangle= F_{1,2}\,  |\varphi_{0} \rangle\,=F_{1,2}\,  \,|1 0 0 \rangle=\frac{1}{\sqrt{3}} \: |1 0 0 \rangle \: + \sqrt{\frac{2}{3}} \: |1 1 0 \rangle $$
NB: Depending on the backend, it happens that the order of the qubits is modified, but without consequence for the state finally reached.
We seem far from the ultimate goal.
Run the following circuit to obtain $|\varphi_{2} \rangle =F_{2,3}\, \, |\varphi_{1} \rangle$
End of explanation
"""


# 3-qubit W state
n = 3
q = QuantumRegister(n) 
c = ClassicalRegister(n)

W_states = QuantumCircuit(q,c) 

W_states.x(q[2]) #start is |100>
F_gate(W_states,q,2,1,3,1) # Applying F12
F_gate(W_states,q,1,0,3,2) # Applying F23

if flag_qx2 : # option ibmqx2 
    W_states.cx(q[1],q[2]) # cNOT 21
    W_states.cx(q[0],q[1]) # cNOT 32
    
else :        # option ibmqx4  
    cxrv(W_states,q,1,2)
    cxrv(W_states,q,0,1)

for i in range(3) :
    W_states.measure(q[i] , c[i]) 

shots = 1024
time_exp = time.strftime('%d/%m/%Y %H:%M:%S')
print('start W state 3-qubit on', backend, "N=", shots,time_exp)
result = execute(W_states, backend=backend, shots=shots)
time_exp = time.strftime('%d/%m/%Y %H:%M:%S')
print('end   W state 3-qubit on', backend, "N=", shots,time_exp)
plot_histogram(result.result().get_counts(W_states))

"""
Explanation: Three-qubit W state, full circuit
In the previous step, we got an histogram compatible with the state:
$$ |\varphi_{2} \rangle =F_{2,3}\,  \, |\varphi_{1} \rangle=F_{2,3}\,  \, (\frac{1}{\sqrt{3}} \: |1 0 0 \rangle \: + \sqrt{\frac{2}{3}} \: |1 1 0 )= \frac{1}{\sqrt{3}} \: (|1 0 0 \rangle \: +  |1 1 0 \:\rangle +  |1 1 1\rangle) $$
NB: Again, depending on the backend, it happens that the order of the qubits is modified, but without consequence for the state finally reached.
It looks like we are nearing the goal.
Indeed, two $cNOT$ gates will make it possible to create a W state.
Run the following cell to see what happens. Did we succeed?
End of explanation
"""


# 4-qubit W state
n = 4
q = QuantumRegister(n) 
c = ClassicalRegister(n)
W_states = QuantumCircuit(q,c) 

W_states.x(q[3]) #start is |1000>
F_gate(W_states,q,3,2,4,1) # Applying F12
F_gate(W_states,q,2,1,4,2) # Applying F23
F_gate(W_states,q,1,0,4,3) # Applying F34        
cxrv(W_states,q,2,3) # cNOT 21
if flag_qx2 : # option ibmqx2  
    W_states.cx(q[1],q[2]) # cNOT 32
    W_states.cx(q[0],q[1]) # cNOT 43

else :        # option ibmqx4 
    cxrv(W_states,q,1,2)
    cxrv(W_states,q,0,1)

for i in range(4) :
    W_states.measure(q[i] , c[i]) 

# circuits = ['W_states']

shots = 1024
time_exp = time.strftime('%d/%m/%Y %H:%M:%S')
print('start W state 4-qubit ', backend, "N=", shots,time_exp)
result = execute(W_states, backend=backend, shots=shots)
time_exp = time.strftime('%d/%m/%Y %H:%M:%S')
print('end   W state 4-qubit on', backend, "N=", shots,time_exp)
plot_histogram(result.result().get_counts(W_states))

"""
Explanation: Now you get an histogram compatible with the final state $|W_{3}\rangle$ through the following steps:
$$ |\varphi_{3} \rangle =  cNOT_{2,1}\,  \, |\varphi_{2} \rangle =cNOT_{2,1}\,\frac{1}{\sqrt{3}} \: (|1 0 0 \rangle \: +  |1 1 0 \rangle\: +  |1 1 1\rangle)  =  \frac{1}{\sqrt{3}} \: (|1 0 0 \rangle \: +  |0 1 0 \: +  |0 1 1\rangle) $$
$$ |W_{3} \rangle =  cNOT_{3,2}\, \, |\varphi_{3} \rangle =cNOT_{3,2}\,\frac{1}{\sqrt{3}} \: (|1 0 0 \rangle \: +  |010 \: \rangle+  |0 1 1\rangle)  =  \frac{1}{\sqrt{3}} \: (|1 0 0 \rangle \: +  |0 1 0 \: +  |0 0 1\rangle) $$
Bingo!
Four-qubit W state
In this section, the production of a four-qubit W state will be obtained by extending the previous circuit.
In this circuit, the starting state is now: $ |\varphi_{0} \rangle \, = \, |1000\rangle$.
A $F$ gate was added at the beginning of the circuit and a $cNOT$ gate was added before the measurement phase.
The entire circuit corresponds to:
$$ |W_{4}\rangle \,=\, cNOT_{4,3}\,  \, cNOT_{3,2}\,  \, cNOT_{2,1}\,  \, F_{3,4} \, \, F_{2,3} \,  \, F_{1,2} \,  \,|\varphi_{0} \rangle \, $$  
Run the following circuit and see what happens.
End of explanation
"""


# 5-qubit W state
n = 5
q = QuantumRegister(n) 
c = ClassicalRegister(n)
W_states = QuantumCircuit(q,c) 

W_states.x(q[4]) #start is |10000>
F_gate(W_states,q,4,3,5,1) # Applying F12
F_gate(W_states,q,3,2,5,2) # Applying F23
F_gate(W_states,q,2,1,5,3) # Applying F34
F_gate(W_states,q,1,0,5,4) # Applying F45
W_states.cx(q[3],q[4]) # cNOT 21
cxrv(W_states,q,2,3)   # cNOT 32
if flag_qx2 : # option ibmqx2 
    W_states.cx(q[1],q[2]) # cNOT 43
    W_states.cx(q[0],q[1]) # cNOT 54
    
else :        # option ibmqx4 
    cxrv(W_states,q,1,2)
    cxrv(W_states,q,0,1)     
    
for i in range(5) :
    W_states.measure(q[i] , c[i]) 

shots = 1024
time_exp = time.strftime('%d/%m/%Y %H:%M:%S')
print('start W state 5-qubit on', backend, "N=", shots,time_exp)
result = execute(W_states, backend=backend, shots=shots)
time_exp = time.strftime('%d/%m/%Y %H:%M:%S')
print('end   W state 5-qubit on', backend, "N=", shots,time_exp)
plot_histogram(result.result().get_counts(W_states))

"""
Explanation: Now, if you used a simulator, you get an histogram clearly compatible with the state:
$$ |W_{4}\rangle \;=\; \frac{1}{2} \: (\:|1000\rangle  +  |0100\rangle  + |0010\rangle + |0001\rangle \:) $$
If you used a real quantum computer, the columns of the histogram compatible with a $|W_{4}\rangle$ state are not all among the highest one. Errors are spreading...
Five-qubit W state
In this section, a five-qubit W state will be obtained, again by extending the previous circuit.
In this circuit, the starting state is now: $ |\varphi_{0} \rangle  =  |10000\rangle$.
A $F$ gate was added at the beginning of the circuit and an additionnal $cNOT$ gate was added before the measurement phase.
$$ |W_{5}\rangle = cNOT_{5,4}   cNOT_{4,3}   cNOT_{3,2}   cNOT_{2,1}  F_{4,5}   F_{3,4}   F_{2,3}    F_{1,2}   |\varphi_{0} \rangle  $$
Run the following cell and see what happens.
End of explanation
"""


