%run init.ipynb

from qiskit import *
nshots = 8192
IBMQ.load_account()
provider= qiskit.IBMQ.get_provider(hub='ibm-q-research-2',group='federal-uni-sant-1',project='main')
device = provider.get_backend('ibmq_lima')
simulator = Aer.get_backend('qasm_simulator')
from qiskit.visualization import plot_histogram
from qiskit.tools.monitor import job_monitor
from qiskit.quantum_info import state_fidelity

def qc_qft(n):
    qr = QuantumRegister(n); qc = QuantumCircuit(qr, name = 'QFT')
    for l in range(0, n):
        qc.h(qr[l])
        if l < n-1:
            for q in range(l+1, n):
                lb = 2*math.pi*2**(-q+l-1)
                qc.cp(lb, qr[q], qr[l])
            #qc.barrier()
    #qc.barrier()
    if n%2 == 0:
        ul = n//2
    else:
        ul = (n-1)//2
    for p in range(0, ul):
        qc.swap(p, n-1-p)
    return qc

n = 4; qc_qft_ = qc_qft(n); qc_qft_.draw(output='mpl')

def qc_iqft(n):
    qr = QuantumRegister(n); qc = QuantumCircuit(qr, name = 'IQFT')
    if n%2 == 0:
        ul = n//2
    else:
        ul = (n-1)//2
    for p in range(ul-1, -1, -1):
        qc.swap(p, n-1-p)
    #qc.barrier()
    for l in range(n-1, -1, -1):
        if l < n-1:
            for q in range(n-1, l+1-1, -1):
                lb = -2*math.pi*2**(-q+l-1)#; print(lb)
                qc.cp(lb, qr[q], qr[l])
        qc.h(qr[l])
        #qc.barrier()
    return qc

n = 4; qc_iqft_ = qc_iqft(n); qc_iqft_.draw(output='mpl')

#x0 = 1 + 0*1j; x1 = 0 + 0*1j # estado |0>
#x0 = 0 + 0*1j; x1 = 1 + 0*1j # estado |1>
x0 = 1/math.sqrt(2) + 0*1j; x1 = 1/math.sqrt(2) + 0*1j # estado |+>

mx0 = math.sqrt(x0.real**2 + x0.imag**2)
if x0.real != 0:
    ph0 = math.atan(x0.imag/x0.real)
elif x0.real == 0 and x0.imag != 0:
    ph0 = math.pi/2
elif x0.real == 0 and x0.imag == 0:
    ph0 = 0
mx1 = math.sqrt(x1.real**2 + x1.imag**2)
if x1.real != 0:
    ph1 = math.atan(x1.imag/x1.real)
elif x1.real == 0 and x1.imag != 0:
    ph1 = math.pi/2
elif x1.real == 0 and x1.imag == 0:
    ph1 = 0
print('|x0|=',mx0,', ph0=', ph0,', |x1|=', mx1,', ph1=', ph1)

th = 2*math.acos(mx1); ph = ph1-ph0+math.pi; lb = ph0-math.pi
print('th=',th,', ph=', ph,', lb=', lb)

n = 1
qr = QuantumRegister(n); qc = QuantumCircuit(qr)
qc.x(qr[0]); qc.u(th, ph, lb, qr[0]); qc.barrier()
qc_qft_ = qc_qft(n); qft = qc_qft_.to_gate(); qc.append(qft, [0])
qc.draw(output='mpl')

svsimulator = Aer.get_backend('statevector_simulator')
job = execute(qc, backend = svsimulator, shots = 1, memory = True)
print(job.result().get_statevector(qc))

def tfc(x): # classical discrete Fourier transform
    d = len(x)
    y = np.zeros(d, dtype = complex)
    for k in range(0, d):
        for j in range(0, d):
            y[k] += x[j]*(math.cos(2*math.pi*j*k/d) + 1j*math.sin(2*math.pi*j*k/d))
    return y/math.sqrt(d)

d = 4; x = np.zeros(d, dtype = complex); 
#x = [1/2, 1/2, 1/2, 1/2]
x = [-1/2, -1/2, 1/2, 1/2]
y = tfc(x); print(y)

n = 2
qr = QuantumRegister(n); qc = QuantumCircuit(qr)
#qc.h([0,1]); qc.barrier() # ok
qc.x([0]); qc.h([0,1]); qc.barrier() # ok a menos de uma fase global de pi
qc_qft_ = qc_qft(n); qft = qc_qft_.to_gate(); qc.append(qft, [0,1])
qc.draw(output='mpl')

svsimulator = Aer.get_backend('statevector_simulator')
job = execute(qc, backend = svsimulator, shots = 1, memory = True)
print(job.result().get_statevector(qc)) # para >1 qubit tem a questao da ordem da bc

d = 4; x = np.zeros(d, dtype = complex); 
x = [math.sqrt(1/6), math.sqrt(1/6), math.sqrt(2/3), 0]
y = tfc(x); print(y)

n = 3
qr = QuantumRegister(n); cr = ClassicalRegister(n); qc = QuantumCircuit(qr, cr)
qc.x([0,2]); qc.barrier()
qc_qft_ = qc_qft(n); qft = qc_qft_.to_gate(); qc.append(qft, [0,1,2]); qc.barrier()
qc_iqft_ = qc_iqft(n); iqft = qc_iqft_.to_gate(); qc.append(iqft, [0,1,2]); qc.barrier()
qc.measure([0,1,2],[0,1,2])
qc.draw(output='mpl')

job_sim = execute(qc, backend = simulator, shots = nshots)
job_exp = execute(qc, backend = device, shots = nshots)
job_monitor(job_exp)

plot_histogram([job_sim.result().get_counts(), job_exp.result().get_counts()], 
               legend = ['sim', 'exp'])

