from qiskit import QuantumCircuit, execute, transpile, Aer 
from qiskit.extensions import UnitaryGate,Initialize
from qiskit.quantum_info import Statevector 
from qiskit.tools.visualization import plot_bloch_vector
from qiskit.tools.visualization import plot_histogram,plot_bloch_multivector  
import numpy as np 
from time import sleep 
import sys 
from scipy.stats import unitary_group 
import matplotlib.pyplot as plt 
%matplotlib inline 
from qiskit import IBMQ 

# defining the absolute error in the eigenvalue
deltas = [0.4, 0.3, 0.2, 0.15, 0.1, 0.05, 0.025, 0.01]

# find the target costs for each of the algorithm 
target_costs = []
for i in deltas:
    cost = (np.cos(i))**2 
    target_costs.append(cost)
target_costs

u1 = unitary_group.rvs(2)
u1 

eigen_phases1, eigen_vectors1 = np.linalg.eig(u1)
print("Eigen states of Unitary :",eigen_vectors1)
eigen_phases1 = np.angle(eigen_phases1)
ep1 = []
for k in eigen_phases1:
    if k < 0:
        ep1.append((k + 2*np.pi)/(2*np.pi))
    else:
        ep1.append(k/(2*np.pi))
eigen_phases1 = np.array(ep1)
print("Eigen phases of unitary :",eigen_phases1)

from normal_SPEA import SPEA 
backend = Aer.get_backend('qasm_simulator')

# looping and storing the results we get 
generated_deltas, errors = [],[]
for delta,cost in zip(deltas,target_costs):
    spea = SPEA(u1,resolution=35,max_iters=15)
    max_differences = []
    delta_error = []
    print("TARGET COST :",cost)
    print("GIVEN DELTA :",delta)
    # doing numerous experiments to converge 
    for exp in range(10):
        result = spea.get_eigen_pair(backend = backend,progress = False, randomize = True, target_cost = cost)
        print("Result with target cost as",cost,"was returned\n",result)
        # find the final cost 
        theta = result['theta']
        stop_cost = result['cost']

        # find the maximum difference between the eigenphases and this returned phase 
        max_diff = max(abs(eigen_phases1 - theta)) 
        max_differences.append(max_diff)

        # now append this and the error from how it actually should have been 
        # NOTE - there is no error if the max difference is actually less than 
        # the delta desired .
        if delta - max_diff >= 0:
            error = 0 
        else:
            # this error should not be more than 5-6 % or 0.05 - 0.06 
            error = abs((delta - max_diff))/delta
        delta_error.append(error)
    # the maximum would be a better metric 
    generated_deltas.append(np.average(max_differences))
    errors.append(np.average(delta_error))


iters = [i for i in range(1,9)]
plt.figure(figsize=(10,8))
plt.title("Plot for Generated and Given maximum $ \delta^{*}$ values",fontsize=16)
plt.xlabel("Iteration Number", fontsize = 14)
plt.ylabel("Values of $ \delta^{*}$ (maximum average error)",fontsize = 14)
plt.grid()
plt.plot(iters, deltas, label = 'Actual Values', color = 'black', alpha = 0.7, marker = 'o')
plt.plot(iters, generated_deltas, label = 'Generated Values', color = 'green', alpha = 0.7, marker = 's')
plt.legend()
plt.savefig("Plot for verifying Cost Lower Bound.JPG",dpi = 200)

plt.figure(figsize=(10,8))
plt.title("Plot for Error in Generated and Given maximum $ \delta^{*}$ values",fontsize=16)
plt.xlabel("Iteration Number", fontsize = 14)
plt.ylabel("Relative Error in $ \delta^{*}$",fontsize = 14)
plt.grid()
plt.plot(iters, errors, label = 'Error Values', color = 'cyan', alpha = 0.8, marker = 'o')
plt.legend()
plt.savefig("Plot for Error in Cost Lower Bound.JPG",dpi = 200)

