import qutip as qt # I use QuTip for the tensor and trace function
import numpy as np 

H = qt.Qobj(np.matrix([[1.0   , 0.0 , 0.0 , 0.0 ]
                      ,[0.0   , 0.0 ,-1.0 , 0.0 ]
                      ,[0.0   ,-1.0 , 0.0 , 0.0 ] 
                      ,[0.0   , 0.0 , 0.0 , 1.0 ]]))

H

Id = qt.qeye(2); Sx = qt.sigmax(); Sy = qt.sigmay(); Sz = qt.sigmaz();

Pauli = [Id,Sx,Sy,Sz]
PauliLabels = ['Id' , 'Sx' , 'Sy' , 'Sz'] 
aij  = 0;

for i in range(4):
    for j in range(4):
        PauliTensor = qt.Qobj(qt.tensor(Pauli[i], Pauli[j]).data.toarray(),dims=[[4],[4]]) # Tensor product of the pauli matrices (sigma_i (tensor) sigma_j), dims[[4],[4]] equal to dimens of  H
        aij = 0.25*( (PauliTensor*H ).tr()) #Elements of the descomposition
        label1 = PauliLabels[i] #Label i to print
        label2 = PauliLabels[j] #Label j to print
        if aij!= 0:
            print(aij, '*', label1, '(tensor)', label2 ) # Print the descomposition of the matrix 

Id_tensor_Id = qt.tensor(Pauli[0], Pauli[0]) 
sx_tensor_sx = qt.tensor(Pauli[1], Pauli[1]) 
sy_tensor_sy = qt.tensor(Pauli[2], Pauli[2]) 
sz_tensor_sz = qt.tensor(Pauli[3], Pauli[3]) 

H_descomposition = 0.5*Id_tensor_Id - 0.5*sx_tensor_sx  - 0.5*sy_tensor_sy + 0.5*sz_tensor_sz
H_descomposition_Dim = qt.Qobj(H_descomposition.data.toarray(),dims=[[4],[4]])

H_descomposition_Dim



from qiskit import QuantumCircuit, Aer, execute 
from qiskit.visualization import plot_histogram
import matplotlib.pyplot as plt

def ansatz(circuit,angle):
    circuit.h(0);         circuit.id(1) #Apply (H \otimes I) |00>
    circuit.cx(0,1)                     #Apply (CX) 
    circuit.rx(angle, 0); circuit.id(1) #Apply (Rx \otimes I)
    return circuit

ansatzCircuit = QuantumCircuit(2)
ansatzCircuit = ansatz(ansatzCircuit,np.pi/4)
ansatzCircuit.draw()

def zz_meas(angle):
    qc = QuantumCircuit(2) # qc is the quantum register of the qubits, in this case 2 qubits
    qc = ansatz(qc,angle)
    qc.measure_all()
    return qc 

zz = zz_meas(np.pi/4)
zz.draw()

def yy_meas(angle):
    qc = QuantumCircuit(2)  # qc is the quantum register of the qubits, in this case 2 qubits
    qc = ansatz(qc,angle)
    qc.barrier()
    
    qc.h(0);  qc.h(1)   #Apply (H \otimes H) |00>
    qc.sdg(0); qc.sdg(1) #Apply (S^{\dagger} \otimes S^{\dagger}) 
    qc.measure_all()
    return qc 

yy = yy_meas(np.pi/4)
yy.draw()

def xx_meas(angle):
    qc = QuantumCircuit(2) # qc is the quantum register of the qubits, in this case 2 qubits
    qc = ansatz(qc,angle)
    
    qc.barrier()

    qc.h(0); qc.h(1) #Apply (H \otimes H) |00>
    qc.measure_all()
    return qc 

xx = xx_meas(np.pi/4)
xx.draw()

def expectation(qcGiven):
    simulator = Aer.get_backend('qasm_simulator') #Simulator for the measures
    result = execute(qcGiven, backend = simulator, shots=10000).result() # Execute the quantum circuit (qcGiven) 10000 times
    counts = result.get_counts(qcGiven) #Count the results
    # 'if' for the times that there are not measure in the state [3]
    if '00' not in counts:
        counts['00'] = 0
    if '01' not in counts:
        counts['01'] = 0
    if '10' not in counts:
        counts['10'] = 0
    if '11' not in counts:
        counts['11'] = 0
       
    total_counts = counts['00'] + counts['11'] + counts['01'] + counts['10']
    zzCounts = counts['00'] + counts['11'] - counts['01'] - counts['10']
    zzCounts = zzCounts / total_counts
    #plot_histogram(counts)
    
    return zzCounts

E = []; # List where the estimations of the eigenvalue will be stored
anglelist = np.linspace(0.0,2*np.pi,200) # Parameter \theta for the minimization

def Eigenvalue_Approx(xx,yy,zz):
    xx_value = expectation(xx) 
    yy_value = expectation(yy)
    zz_value = expectation(zz)
    EigenValue = 0.5 - (0.5)*xx_value - (0.5)*yy_value + (0.5)*zz_value # Estimate the expectation value of the matriz H
    return EigenValue 

for theta in anglelist:
    xx = xx_meas(theta)
    yy = yy_meas(theta)
    zz = zz_meas(theta)
    E.append(Eigenvalue_Approx(xx,yy,zz))

plt.figure(figsize = (10,4))
plt.plot(anglelist,E , linewidth=5)
plt.grid()


plt.xlabel(r'$\theta$', fontsize = 22)
plt.ylabel(r'$<{\psi(\theta)} | H |{\psi(\theta)} >$' , fontsize = 20)
plt.yticks(fontsize = 16)
plt.xticks([0.0, np.pi/4, np.pi/2, 3*np.pi/4 , np.pi, 5*np.pi/4 , 3*np.pi/2, 7*np.pi/4 , 2*np.pi],
           [r'$0$', r'$\pi/4$', r'$\pi/2$', r'$3\pi/4$', r'$\pi$', r'$5\pi/4$', r'$3\pi/2$' , r'$7\pi/4$', r'$2\pi$' ], fontsize = 16)
plt.show()

MimimumEigenvalue = np.amin(E)
print('The estimated minimum eigenvalue of the matrix H from VQE algorithm is: {}'.format(MimimumEigenvalue))

evalsH, eketsH =H.eigenstates()
MimimumEigenvalue_Exact  = np.amin(evalsH)
print('The exact minimum eigenvalue of the matrix H is: {}'.format(MimimumEigenvalue_Exact))

evalsH

evalsH[0]

