from sympy import *
init_printing(use_unicode=True)
%matplotlib inline

p00,p01,p10,p11 = symbols('p_{00} p_{01} p_{10} p_{11}')
th,ph = symbols('theta phi')
Psi00 = Matrix([[cos(th)],[0],[0],[sin(th)]])
Psi01 = Matrix([[0],[sin(ph)],[cos(ph)],[0]])
Psi11 = Matrix([[0],[cos(ph)],[-sin(ph)],[0]])
Psi10 = Matrix([[-sin(th)],[0],[0],[cos(th)]])
#Psi00, Psi00.T, Psi00*Psi00.T

rhoX = p00*Psi00*Psi00.T +  p01*Psi01*Psi01.T + p10*Psi10*Psi10.T + p11*Psi11*Psi11.T
simplify(rhoX)

def kp(x,y):
    return KroneckerProduct(x,y)
I = Matrix([[1,0],[0,1]])
Y = Matrix([[0,-1j],[1j,0]])
Y = Matrix([[0,1],[1,0]])
Z = Matrix([[1,0],[0,-1]])
cxx,cyy,czz,az,bz = symbols('c_{xx} c_{yy} c_{zz} a_{z} b_{z}')
rhoX = (1/4)*(kp(I,I) + cxx*kp(X,X) + cyy*kp(Y,Y) + czz*kp(Z,Z) + az*kp(Z,I) + bz*kp(I,Z))
simplify(rhoX)

th,be,ga = symbols('theta beta gamma')
c00 = cos(th/2)*cos(be/2)*cos(ga/2) + sin(th/2)*sin(be/2)*sin(ga/2) 
c01 = cos(th/2)*cos(be/2)*sin(ga/2) - sin(th/2)*sin(be/2)*cos(ga/2)
c10 = cos(th/2)*sin(be/2)*cos(ga/2) - sin(th/2)*cos(be/2)*sin(ga/2)
c11 = cos(th/2)*sin(be/2)*sin(ga/2) + sin(th/2)*cos(be/2)*cos(ga/2)

simplify(c00**2 + c01**2 + c10**2 + c11**2) # ok!

P0 = Matrix([[1,0],[0,0]])
P1 = Matrix([[0,0],[0,1]])
def Ry(th):
    return cos(th/2)*I - 1j*sin(th/2)*Y
def Cx_ab():
    return KroneckerProduct(P0,I) + KroneckerProduct(P1,X) 
def Cx_ba():
    return KroneckerProduct(I,P0) + KroneckerProduct(X,P1) 

MB = Cx_ab()*KroneckerProduct(Ry(th-ph),I)*Cx_ba()*KroneckerProduct(Ry(th+ph),I) # mudanca de base
simplify(MB)

from qiskit import *
import numpy as np
import math
import qiskit
nshots = 8192
IBMQ.load_account()
provider = qiskit.IBMQ.get_provider(hub='ibm-q', group='open', project='main')
device = provider.get_backend('ibmq_belem')
simulator = Aer.get_backend('qasm_simulator')
from qiskit.tools.monitor import job_monitor
from qiskit.ignis.verification.tomography import state_tomography_circuits, StateTomographyFitter
from qiskit.ignis.mitigation.measurement import complete_meas_cal, CompleteMeasFitter

# retorna o circuito quantico que prepara um certo estado real de 1 qubit
# coef = array com os 2 coeficientes reais do estado na base computacional
def qc_psi_1qb_real(coef):
    qr = QuantumRegister(1)
    qc = QuantumCircuit(qr, name='psir_1qb')
    th = 2*math.acos(np.abs(coef[0]))
    qc.ry(th, qr[0])
    return qc

eigvals = [0.1,0.9]
coef = np.sqrt(eigvals)
print(coef)
qc_psi_1qb_real_ = qc_psi_1qb_real(coef)
qc_psi_1qb_real_.draw('mpl')

from qiskit.quantum_info import Statevector
sv = Statevector.from_label('0')
sv

sv = sv.evolve(qc_psi_1qb_real_)
sv

# retorna o circuito quantico que prepara um certo estado real de 2 qubits
# coef = array com os 4 coeficientes reais do estado na base computacional
def qc_psi_2qb_real(coef):
    qr = QuantumRegister(2)
    qc = QuantumCircuit(qr, name = 'psir_2qb')
    xi = 2*math.acos(math.sqrt(coef[0]**2+coef[1]**2))
    coef1 = [math.sqrt(coef[0]**2+coef[1]**2),math.sqrt(1-coef[0]**2-coef[1]**2)]
    c_psi_1qb_real_ = qc_psi_1qb_real(coef1)
    qc.append(c_psi_1qb_real_, [qr[0]])
    th0 = 2*math.atan(np.abs(coef[1])/np.abs(coef[0]))
    th1 = 2*math.atan(np.abs(coef[3])/np.abs(coef[2]))
    qc.x(0)
    qc.cry(th0, 0, 1)
    qc.x(0)
    qc.cry(th1, 0, 1)
    return qc

eigvals = [0.1, 0.2, 0.3, 0.4]
coef = np.sqrt(eigvals)
print(coef)
qc_psi_2qb_real_ = qc_psi_2qb_real(coef)
qc_psi_2qb_real_.draw('mpl')

sv = Statevector.from_label('00')
sv

sv = sv.evolve(qc_psi_2qb_real_)
sv # note o ordenamento, o 2º qb é o 1º e o 1º é o 2º (00,10,01,11)

def qc_ry(th):
    qr = QuantumRegister(1)
    qc = QuantumCircuit(qr, name = 'RY')
    qc.ry(th, 0)
    return qc

# retorna o circuito quantico que prepara um certo estado real de 3 qubits
# coef = array com os 8 coeficientes reais do estado na base computacional
def qc_psi_3qb_real(coef):
    qr = QuantumRegister(3)
    qc = QuantumCircuit(qr, name = 'psir_3qb')
    d = len(coef)
    coef2 = np.zeros(d//2)
    th = np.zeros(d//2)
    for j in range(0,2):
        for k in range(0,2):
            coef2[int(str(j)+str(k),2)] = math.sqrt(coef[int(str(j)+str(k)+str(0),2)]**2+coef[int(str(j)+str(k)+str(1),2)]**2) 
    c_psi_2qb_real_ = qc_psi_2qb_real(coef2)
    qc.append(c_psi_2qb_real_, [qr[0],qr[1]])
    for j in range(0,2):
        for k in range(0,2):
            th[int(str(j)+str(k),2)] = 2*math.atan(np.abs(coef[int(str(j)+str(k)+str(1),2)])/np.abs(coef[int(str(j)+str(k)+str(0),2)]))
            if j == 0:
                qc.x(0)
            if k == 0:
                qc.x(1)
            qc_ry_ = qc_ry(th[int(str(j)+str(k),2)])
            ccry = qc_ry_.to_gate().control(2)
            qc.append(ccry, [0,1,2])
            if j == 0:
                qc.x(0)
            if k == 0:
                qc.x(1)
    return qc

list_bin = []
for j in range(0,2**3):
    b = "{:03b}".format(j)
    list_bin.append(b)
print(list_bin)

eigvals = [0.01,0.1,0.04,0.2,0.05,0.3,0.18,0.12]
coef = np.sqrt(eigvals)
print(coef)
qc_psi_3qb_real_ = qc_psi_3qb_real(coef)
qc_psi_3qb_real_.draw('mpl') # odernamento: '000', '001', '010', '011', '100', '101', '110', '111'

sv = Statevector.from_label('000')
sv

sv = sv.evolve(qc_psi_3qb_real_)
sv # ordenamento aqui: 000 100 010 110 001 101 011 111

# retorna o circuito quantico que prepara um certo estado real de 4 qubits
# coef = array com os 16 coeficientes reais do estado na base computacional
def qc_psi_4qb_real(coef):
    qr = QuantumRegister(4)
    qc = QuantumCircuit(qr, name = 'psir_4qb')
    d = len(coef)
    coef3 = np.zeros(d//2)
    th = np.zeros(d//2)
    for j in range(0,2):
        for k in range(0,2):
            for l in range(0,2):
                coef3[int(str(j)+str(k)+str(l),2)] = math.sqrt(coef[int(str(j)+str(k)+str(l)+str(0),2)]**2+coef[int(str(j)+str(k)+str(l)+str(1),2)]**2) 
    c_psi_3qb_real_ = qc_psi_3qb_real(coef3)
    qc.append(c_psi_3qb_real_, [qr[0],qr[1],qr[2]])
    for j in range(0,2):
        for k in range(0,2):
            for l in range(0,2):
                th[int(str(j)+str(k)+str(l),2)] = 2*math.atan(np.abs(coef[int(str(j)+str(k)+str(l)+str(1),2)])/np.abs(coef[int(str(j)+str(k)+str(l)+str(0),2)]))
                if j == 0:
                    qc.x(0)
                if k == 0:
                    qc.x(1)
                if l == 0:
                    qc.x(2)
                qc_ry_ = qc_ry(th[int(str(j)+str(k)+str(l),2)])
                ccry = qc_ry_.to_gate().control(3)
                qc.append(ccry, [0,1,2,3])
                if j == 0:
                    qc.x(0)
                if k == 0:
                    qc.x(1)                
                if l == 0:
                    qc.x(2)
    return qc

list_bin = []
for j in range(0,2**4):
    b = "{:04b}".format(j)
    list_bin.append(b)
print(list_bin)

eigvals = np.zeros(2**4)
eigvals[0] = 0.008
for j in range(1,len(eigvals)-1):
    eigvals[j] = eigvals[j-1]+0.005
#print(np.sum(eigvals))
eigvals[j+1] = 1 - np.sum(eigvals)
#print(eigvals)
#print(np.sum(eigvals))
coef = np.sqrt(eigvals)
print(coef)
qc_psi_4qb_real_ = qc_psi_4qb_real(coef)
qc_psi_4qb_real_.draw('mpl') 
# '0000', '0001', '0010', '0011', '0100', '0101', '0110', '0111', '1000', '1001', '1010', '1011', '1100', '1101', '1110', '1111'

sv = Statevector.from_label('0000')
sv

sv = sv.evolve(qc_psi_4qb_real_)
sv # ordenamento aqui: 0000 1000 0100 1100 0010 1010 0110 1110 0001 1001 0101 1101 0011 1011 0111 1111

sv[1]









