%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 executing 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 executing V
qc2.x(0)
qc2.toffoli(0,1,2)#Flipping the third bit if the first is zero and second is one
qc2.x(0)          #Essentially 010 -> 011

qc2.toffoli(1,2,0)

qc2.x(0)
qc2.toffoli(0,1,2)#Undoing the flip from before
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)
qcy.t(0)
qcy.t(0)
qcy.t(0)
qcy.t(0)
qcy.h(0)
qcy.t(0)
qcy.t(0)
qcy.t(0)
qcy.t(0)
qcy.h(0)# Essentially XZ is being applied since T^4 = X and HXH = Z and XZ = iY but global phase is ignored
Y_circ = qi.Operator(qcy).data
Y = np.array([[0,-1j],[1j,0]])
print(norm(Y_circ,Y))

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

qc3 = QuantumCircuit(1)
qc_min = 1000
while(True):
    parameters = np.random.randint(8, size=6) #Samling integers from 0 to 7 and getting 6 of them
    for i in parameters[0:5]:
        for j in range(i):
            qc3.t(0)
        qc3.h(0)
    for j in range(parameters[-1]):
        qc3.t(0)
    if norm(qi.Operator(qc3).data,uni_q) <= 1e-8:#Break if we have reached the right combination
        print(parameters)
        break
    qc_min = np.sum(parameters)
    qc3 = QuantumCircuit(1)#Reinitiallize circuit

uni_q_circ = qi.Operator(qc3).data
print(norm(uni_q_circ,uni_q))

qc3.draw('mpl')#Drawing the actual solution
#Note that there are multiple solutions so as long as your norm
#was of order 1.65e-9 you are correct. The reason of these multiple
#solutions is because there are commutativity relations that we
#havent bothered calculating

#There are better ways and many of you probably did use them. This is the best circuit we found
#If you search amongst all the possibilities or iterate the right way this will be the solution you find
#The difference between this cell and the one above is that the one above came from a random number generator
qc3 = QuantumCircuit(1)
qc3.t(0)
qc3.h(0)
qc3.t(0)
qc3.t(0)
qc3.t(0)
qc3.h(0)
qc3.t(0)
qc3.h(0)
qc3.t(0)
qc3.h(0)
qc3.t(0)
qc3.t(0)
qc3.h(0)
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)

alpha = pi_alpha/np.pi
i = 0
while(True):
    if abs((alpha*i)%2 - 0.5) <= 0.001:
        print(i)
        break
    i += 1

i = 1
while(True):
    if abs((alpha*i)%2) <= 0.013 and (i%4 == 2):
        print(i)
        break
    i += 1

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

#Doing the Rx(pi/2)
for i in range(476):
    qcb.append(controlled_irx,[0,1,2])

#Controlled not off first qubit (here first qubit is actually the third one cause we have two ancillas)
qcb.cx(2,3)

#Doing the Rx(3pi/2) which is Rx(-pi/2) essentially
qcb.x(2)
for i in range(476):
    qcb.append(controlled_irx,[0,1,2])

# This makes the state |11> flip sign if third and fourth qubit is 1 
#so that way we know that fourth qubit has |+> state always
qcb.x(0)#reset ancillas
qcb.x(1)#reset ancillas
qcb.cx(2,0)
qcb.cx(3,1)
for i in range(122):
    qcb.append(controlled_irx,[0,1,3])
qcb.cx(3,1)
qcb.cx(2,0)


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