from qiskit import IBMQ, Aer
from qiskit.aqua import QuantumInstance
from qiskit.providers.aer.noise import NoiseModel
from qiskit.ignis.mitigation.measurement import CompleteMeasFitter

#Loading data about our backend quantum architecture
IBMQ.load_account()
provider = IBMQ.get_provider(hub='ibm-q')
backend_sim = Aer.get_backend("qasm_simulator")
backend_real = provider.get_backend("ibmq_lima")
coupling_map = backend_real.configuration().coupling_map
# Extracting noise model for error mitigation
noise_model = NoiseModel.from_backend(backend_real)
# Initializing quantum instance with error mitigation based on noise model of backend_real
quantum_instance = QuantumInstance(backend=backend_sim,
                                   shots=20000,
                                   noise_model=noise_model,
                                   coupling_map=coupling_map,
                                   measurement_error_mitigation_cls=CompleteMeasFitter,
                                   cals_matrix_refresh_period=30)


from qiskit.aqua.components.optimizers import COBYLA

optimizer = COBYLA(maxiter=5000, tol=0.0001)

from qiskit.aqua.operators import MatrixOperator, op_converter
import numpy as np

#Initialization of a 2x2 matrix in the same shape as MD program
dimensions = 2
rdist = np.random.rand(dimensions,dimensions)
rdist[0][0] = 1.0
rdist[0][1] = 2.0
rdist[1][0] = 1.0
rdist[1][1] = 3.0
rdistt = np.transpose(rdist)
z1 = np.zeros((dimensions,dimensions))
z2 = np.zeros((dimensions,dimensions))
in_matrix = np.block([[z1,rdist],[rdistt,z2]])

#Converting matrix into operator
hamiltonian = MatrixOperator(in_matrix)
hamiltonian_qubit_op = op_converter.to_weighted_pauli_operator(hamiltonian)

from qiskit.circuit.library import RealAmplitudes
best_ansatz_ever = RealAmplitudes(2,reps=2)
best_ansatz_ever.draw(output='text')

from qiskit.aqua.algorithms import VQE

vqe = VQE(operator=hamiltonian_qubit_op,var_form=best_ansatz_ever,quantum_instance=quantum_instance, optimizer=optimizer)
vqe_result = np.real(vqe.run(backend_sim)['eigenvalue'])
print(vqe_result)

import math
import scipy.linalg as lin_alg

#storing values of, respectively, classic and quantum eigenvalues
classic = []
quantum = []

for i in range(1,100):
    #Initializing random matrix
    rdist = np.random.rand(dimensions,dimensions)
    rdistt = np.transpose(rdist)
    z1 = np.zeros((dimensions,dimensions))
    z2 = np.zeros((dimensions,dimensions))
    in_matrix = np.block([[z1,rdist],[rdistt,z2]])
    
    #Conversion of matrix into operator
    hamiltonian = MatrixOperator(in_matrix)
    hamiltonian_qubit_op = op_converter.to_weighted_pauli_operator(hamiltonian)
    
    #Calculation of Quantum eigenvalue
    vqe = VQE(operator=hamiltonian_qubit_op,var_form=best_ansatz_ever,quantum_instance=quantum_instance, optimizer=optimizer)
    vqe_result = np.real(vqe.run(backend_sim)['eigenvalue'])
    
    #Calculation of Classic eigenvalue
    classic_result = lin_alg.eigvalsh(in_matrix)[0]
    
    #Storing eigenvalues of current matrix
    classic.append(classic_result)
    quantum.append(vqe_result)

#Calculating mean square error, root mean square error, root mean square error
MSE = np.square(np.subtract(classic,quantum)).mean()
RMSE = math.sqrt(MSE)
print("RMSE: "+str(RMSE))
format_rmse = "{:.4f}".format((RMSE/(max(classic)-min(classic)))*100.0)
print("NRMSE: "+format_rmse+ "%")

