# Importing standard Qiskit libraries and configuring account
from qiskit import *
from qiskit.compiler import *
from qiskit.tools.jupyter import *
from qiskit.visualization import *
import matplotlib.pyplot as plotter
import scipy
import numpy as np
from IPython.display import display, Math, Latex
import qiskit.quantum_info as qi
%matplotlib inline
%pip install -I git+https://github.com/mnp-club/MnP_QC_Workshop.git
from mnp_qc_workshop_2020.unitary_circuit import *

U = get_unitary()
print(U.shape)
fig, (ax1, ax2) = plotter.subplots(nrows = 1, ncols = 2, figsize=(12,6))
ax1.imshow(np.real(U)) #plot real parts of each element
ax2.imshow(np.imag(U)) #plot imaginary parts of each element
fig, (ax3, ax4) = plotter.subplots(nrows = 1, ncols = 2, figsize=(12,6))
ax3.imshow(np.abs(U)) #plot the absolute values of each element
ax4.imshow(np.angle(U)) #plot the phase angles of each element

qc1 = QuantumCircuit(4)
qc1.unitary(U,range(4))
qc1 = transpile(qc1,basis_gates=['cx','u3'],optimization_level=3)
qc1.draw('mpl')

#Run this cell for getting your circuit checked
check_circuit(qc1)

H = scipy.linalg.hadamard(16)/4 #Normalizing has to be done
HU = np.matmul(H,U)
fig, (ax1, ax2) = plotter.subplots(nrows = 1, ncols = 2, figsize=(12,6))
ax1.imshow(np.real(HU)) #plot real parts of each element
ax2.imshow(np.imag(HU)) #plot imaginary parts of each element
fig, (ax3, ax4) = plotter.subplots(nrows = 1, ncols = 2, figsize=(12,6))
ax3.imshow(np.abs(HU)) #plot the absolute values of each element
ax4.imshow(np.angle(HU)) #plot the phase angles of each element

qc2 = QuantumCircuit(4)
#for i in range(4):
 #   qc2.h(i)
qc2.unitary(HU,range(4))
for i in range(4):
    qc2.h(i)
qc2 = transpile(qc2,basis_gates=['cx','u3'],optimization_level=3)
qc2.draw('mpl')

#Run this cell for getting your circuit checked
check_circuit(qc2)

HUH = np.matmul(HU,H)
fig, (ax1, ax2) = plotter.subplots(nrows = 1, ncols = 2, figsize=(12,6))
ax1.imshow(np.real(HUH)) #plot real parts of each element
ax2.imshow(np.imag(HUH)) #plot imaginary parts of each element
fig, (ax3, ax4) = plotter.subplots(nrows = 1, ncols = 2, figsize=(12,6))
ax3.imshow(np.abs(HUH)) #plot the absolute values of each element
ax4.imshow(np.angle(HUH)) #plot the phase angles of each element

qc3 = QuantumCircuit(4)
for i in range(4):
    qc3.h(i)
qc3.unitary(HUH,range(4))
for i in range(4):
    qc3.h(i)
qc3 = transpile(qc3,basis_gates=['cx','u3'],optimization_level=3)
qc3.draw('mpl')

#Run this cell for getting your circuit checked
check_circuit(qc3)

diag = np.diag(HUH)
print(diag)

a = np.sqrt(0.5)
clean_diag = [ 1, a+(a*1j), -a+(a*1j), -1, 1j, a-(a*1j), -a-(a*1j), 1j, 
                       -a-(a*1j), 1j, -1, a+(a*1j), a-(a*1j), 1, -1j, a-(a*1j)]
print(clean_diag)

qc31 = QuantumCircuit(4)
for i in range(4):
    qc31.h(i)
qc31.diagonal(clean_diag,qc31.qubits)
for i in range(4):
    qc31.h(i)
qc31 = transpile(qc31,basis_gates=['cx','u3'],optimization_level=3)
qc31.draw('mpl')

#Run this cell for getting your circuit checked
check_circuit(qc31)

func = lambda a: np.array([[1,0],[0,np.exp(1j*np.pi*a/4)]])

def checking_of_diag(diag_v):
    for i in range(16):
        if np.abs(diag_v[i]-clean_diag[i])>=0.001 and np.abs(diag_v[i]+clean_diag[i])>=0.001:
            return False
    return True
            

while(True):
    par = np.random.randint(low=1,high=8, size=4)#We will assume that all atleast rotate by pi/4
    V_test = np.kron(np.kron(func(par[0]),func(par[1])),np.kron(func(par[2]),func(par[3])))
    if checking_of_diag(np.diag(V_test)):
        print(par)
        break
    

differences = np.identity(16)
for i in range(16):
    if np.abs(V_test[i,i]-clean_diag[i]) >= 0.0001:
        print(i)
        differences[i,i] = -1

def without_global_phase(matrix: np.ndarray, atol: float = 1e-8) :
    phases1 = np.angle(matrix[abs(matrix) > atol].ravel(order='F'))
    if len(phases1) > 0:
        matrix = np.exp(-1j * phases1[0]) * matrix
    return matrix
#Function just returns norm ignoring global phase between unitaries
def norm(unitary_a: np.ndarray, unitary_b: np.ndarray) :
    return np.linalg.norm(without_global_phase(unitary_b)-without_global_phase(unitary_a), ord=2)



qc_d = QuantumCircuit(4)
for i in range(64):
    if (i%2) == 1:
        qc_d.cz(0,1)
    if (i//2)%2 == 1:
        qc_d.cz(0,2)
    if (i//4)%2 == 1:
        qc_d.cz(0,3)
    if (i//8)%2 == 1:
        qc_d.cz(1,2)
    if (i//16)%2 == 1:
        qc_d.cz(1,3)
    if (i//32)%2 == 1:
        qc_d.cz(3,2)
    if norm(qi.Operator(qc_d).data,differences) <= 0.001:
        print(i)
        break
    qc_d = QuantumCircuit(4)

qc = QuantumCircuit(4)
theta1 = par[3]*np.pi/4
theta2 = par[2]*np.pi/4
theta3 = par[1]*np.pi/4
theta4 = par[0]*np.pi/4
qc.h(0)
qc.h(1)
qc.h(2)
qc.h(3)
qc.rz(theta1,0)
qc.rz(theta2,1)
qc.rz(theta3,2)
qc.rz(theta4,3)
qc.cz(1,3)
qc.cz(0,2)
qc.cz(0,3)
qc.h(0)
qc.h(1)
qc.h(2)
qc.h(3)
qc = transpile(qc,basis_gates=['cx','u3'],optimization_level=3)
qc.draw('mpl')

#Run this cell for getting your circuit checked
check_circuit(qc)