%matplotlib inline
# Importing standard Qiskit libraries and configuring account
from qiskit import *
from qiskit.compiler import *
from qiskit.tools.jupyter import *
from qiskit.visualization import *
import numpy as np
import qiskit.quantum_info as qi
# Loading your IBM Q account(s)
provider = IBMQ.load_account()

#This cell is just for seeing what the U tilde matrix looks like
a, b, c, d = np.cos(np.pi/8), np.sin(np.pi/8), -np.sin(np.pi/8), np.cos(np.pi/8)
u_tilde = np.array([[a,c],[b,d]]).reshape(2,2)
print(u_tilde)

#This cell defines the controlled variant of the U tilde matrix
qc_for_u = QuantumCircuit(1)
qc_for_u.ry(np.pi/4, 0)
qc_for_u.name = "U"
controlled_u_tilde = qc_for_u.to_gate().control(2)

qc1 = QuantumCircuit(3)

qc1.x(0)
qc1.x(1)
qc1.toffoli(0,1,2)#Flipping the third bit if the first two are zero
qc1.x(0)          #Essentially 000 -> 001
qc1.x(1)

qc1.x(0)
qc1.toffoli(0,2,1)#Flipping the second bit if the first bit is zero and the third is one
qc1.x(0)          #Essentially 001 -> 011

qc1.append(controlled_u_tilde, [1, 2, 0])

qc1.x(0)
qc1.toffoli(0,2,1)#Undoing the flip from before
qc1.x(0)

qc1.x(0)
qc1.x(1)
qc1.toffoli(0,1,2)#Undoing the flip from before
qc1.x(0)
qc1.x(1)

qc1.draw('mpl')


U_circ = qi.Operator(qc1).data
print(U_circ)

qc2 = QuantumCircuit(3)

# Code for U

qc2.x(0)
qc2.x(1)
qc2.toffoli(0,1,2)#Flipping the third bit if the first two are zero
qc2.x(0)          #Essentially 000 -> 001
qc2.x(1)

qc2.x(0)
qc2.toffoli(0,2,1)#Flipping the second bit if the first bit is zero and the third is one
qc2.x(0)          #Essentially 001 -> 011

qc2.append(controlled_u_tilde, [1, 2, 0])

qc2.x(0)
qc2.toffoli(0,2,1)#Undoing the flip from before
qc2.x(0)

qc2.x(0)
qc2.x(1)
qc2.toffoli(0,1,2)#Undoing the flip from before
qc2.x(0)
qc2.x(1)

# Code for V

qc2.x(0)
qc2.toffoli(0,1,2) # |010> -> |011>
qc2.x(0)

qc2.toffoli(1,2,0) # V-tilde is the same as Toffoli

qc2.x(0)
qc2.toffoli(0,1,2) # Reversing the flip
qc2.x(0)

qc2.draw('mpl')

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
V_circ = without_global_phase(qi.Operator(qc2).data)
print(V_circ)

#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)

#Make the Pauli Y gate
qcy = QuantumCircuit(1)

for i in range(6):
    qcy.t(0)
    
qcy.h(0)

for i in range(4):
    qcy.t(0)
    
qcy.h(0)

for i in range(2):
    qcy.t(0)

Y_circ = qi.Operator(qcy).data
Y = np.array([[0,-1j],[1j,0]])
print(norm(Y_circ,Y))

print("Final Circuit:")
qcy.draw('mpl')

uni_q = np.array([[-0.25+0.60355339j, 0.60355339+0.45710678j],
                   [0.60355339-0.45710678j, 0.25+0.60355339j]]).reshape(2,2)

while True:
    qc3 = QuantumCircuit(1)
    for i in range(5):
        for j in range(np.random.randint(8)):
            qc3.t(0)
        qc3.h(0)
    for j in range(np.random.randint(8)):
        qc3.t(0)
    uni_q_circ = qi.Operator(qc3).data
    if norm(uni_q_circ,uni_q) < 1e-8:
        break
    
print("Final Error: ", norm(uni_q_circ,uni_q))
print("Final Circuit:")
qc3.draw('mpl')

#Defining the gate as controlled_irx
from qiskit.extensions import *
pi_alpha = np.arccos(0.6)
qc_for_irx = QuantumCircuit(1)
irx = np.array([[1j*np.cos(pi_alpha/2),np.sin(pi_alpha/2)],[np.sin(pi_alpha/2),1j*np.cos(pi_alpha/2)]]).reshape(2,2)
g_irx = UnitaryGate(data=irx,label=None)
controlled_irx = g_irx.control(2)

def final_vector(qc):
    backend = Aer.get_backend('statevector_simulator')
    job = execute(qc, backend)
    result = job.result()
    outputstate = result.get_statevector(qc, decimals=3)
    return outputstate

n = 4
qcb = QuantumCircuit(n,n)
qcb.x(0)#Ancilla
qcb.x(1)#Ancilla  You can use these two as control to use the controlled_irx gate

count1 = 0
while abs(0.5 - (count1*pi_alpha/np.pi)%2) >= 1e-2 or count1 % 4 != 0:
    count1 += 1

count2 = 0
while abs(1.5 - (count2*pi_alpha/np.pi)%2) >= 1e-2 or count2 % 4 != 0:
    count2 += 1

count3 = 0
while abs(4.0 - (count3*pi_alpha/np.pi)%4) >= 1e-2 or count3 % 4 != 2:
    count3 += 1

for i in range(count1):    # Rotation about x by pi/2 (count1 % 4 == 0)
    qcb.append(controlled_irx, [0, 1, 2])
qcb.cx(2, 3)
for i in range(count2):    # Rotation about x by -pi/2 (count2 % 4 == 0)
    qcb.append(controlled_irx, [0, 1, 2])
qcb.x(0)
qcb.cx(3, 0)
for i in range(count3):    # Rotation about x by 0 (count3 % 4 == 2)
    qcb.append(controlled_irx, [0, 2, 3])
qcb.cx(3, 0)
qcb.x(0)

print("Executing circuit...")
        
#Get state of qubit which should have the |+> state using the backend simulator
i = 3 #Index for the qubit at |+> state
qcb.h(i)#Puts the |+> state to |0>

for i in range(4):
    qcb.measure(i, i)

def run_circuit(qcb):
    backend = Aer.get_backend('qasm_simulator') # we choose the simulator as our backend
    result = execute(qcb, backend, shots = 2000).result() # we run the simulation
    counts = result.get_counts() # we get the counts
    return counts

counts = run_circuit(qcb)
print(counts)
plot_histogram(counts)