import numpy as np
import matplotlib.pyplot as plt

# Importing standard Qiskit libraries
from qiskit import QuantumCircuit, transpile, Aer, IBMQ
from qiskit.tools.jupyter import *
from qiskit.visualization import *
# from ibm_quantum_widgetsets import *
from qiskit.providers.aer import QasmSimulator
from qiskit.algorithms.optimizers import *

from qiskit.utils import algorithm_globals
from qiskit.circuit.library import TwoLocal
from qiskit_optimization.algorithms import MinimumEigenOptimizer
from qiskit.algorithms import VQE, QAOA, NumPyMinimumEigensolver
from qiskit.utils import QuantumInstance
from qiskit_optimization import QuadraticProgram

from qiskit.opflow import Z, I
from ibm_quantum_widgets import *
from qiskit.visualization import plot_histogram
# from qiskit.quantum_info.operators import Operator, Pauli

from qiskit.opflow import PauliExpectation, CVaRExpectation
from qiskit.circuit.library import RealAmplitudes
from qiskit_optimization.converters import LinearEqualityToPenalty

from qiskit_optimization.algorithms import GroverOptimizer, MinimumEigenOptimizer
# Loading your IBM Quantum account(s)
provider = IBMQ.load_account()

def output(result, prob_path, fval_path):
    prob = {}
    fval = {}
    for i in range(len(result.samples)):
        tmp = ""
        count = 0
        for num in range(8):
            tmp += str(int(result.samples[i].x[num]))
            count += int(result.samples[i].x[num])

        prob[tmp] = result.samples[i].probability
        fval[tmp] = result.samples[i].fval


    np.save(prob_path, prob)
    np.save(fval_path, fval)

def unnormalized_h(adj_mtx):
    qubit_num = len(adj_mtx)
    
    iden = I
    for i in range(1,qubit_num):
        iden = iden^I
    
    op = iden-iden
    # print(op)
    
    for i in range(qubit_num):
        for j in range(qubit_num):
            if i > j:
                # 2ZiZj
                temp = np.ones(qubit_num)*I
                temp[i] = Z
                temp[j] = Z
                
                op_0 = temp[0]
                for k in range(1, qubit_num):
                    op_0 = op_0^temp[k]
                    
                print(op_0)
                op = op + 2*op_0
                
                # 0.5*I - 0.5*ZiZj
                if adj_mtx[i][j] == 1:
                    op = op + (0.5*iden - 0.5 * op_0)
    op = op + qubit_num * iden
        
    return op

def normalized_h(adj_mtx,ratio):
    qubit_num = len(adj_mtx)
    
    iden = I
    for i in range(1,qubit_num):
        iden = iden^I
    
    op = iden-iden
    #print(op)
    
    # for 1/m * 0.5*sum(I-zizj) term
    first_term = iden-iden
    m=0
    
    # for zizj term
    second_term = iden - iden
    
    for i in range(qubit_num):
        for j in range(qubit_num):
            if i>j:
                # 2ZiZj
                temp = np.ones(qubit_num)*I
                temp[i] = Z
                temp[j] = Z
                
                op_0 = temp[0]
                for k in range(1,qubit_num):
                    op_0 = op_0^temp[k]
                    
                #print(op_0)
                second_term = second_term + 2*op_0
                
                # 0.5*I - 0.5*ZiZj
                if adj_mtx[i][j] == 1:
                    first_term = first_term + 0.5*iden - 0.5 * op_0
                    m+=1
    
    #unweighted_h = op + first_term + second_term + qubit_num * iden
    k = ((ratio/qubit_num**2))
    weighted_h = op + first_term/m + float(k)*(second_term + qubit_num * iden)
        
    return weighted_h

def adjacency_matrix(graph):
    matrix = []
    for i in range(len(graph)): 
        matrix.append([0]*(len(graph)))
        for j in graph[i]:
            matrix[i][j] = 1
    return matrix

lst = [[1,4,6],[0,2],[1,5,7],[4],[0,3,5],[2,4],[0],[2]]
adj_mtx = adjacency_matrix(lst)
op = unnormalized_h(adj_mtx)
# op.to_matrix_op()  # See the operator with matrix form
qp = QuadraticProgram()
qp.from_ising(op)
print(qp)

# solving Quadratic Program using exact classical eigensolver
exact = MinimumEigenOptimizer(NumPyMinimumEigensolver())
result = exact.solve(qp)
print(result.prettyprint())

lst = [[1,4,6],[0,2],[1,5,7],[4],[0,3,5],[2,4],[0],[2]]
adj_mtx = adjacency_matrix(lst)
op = unnormalized_h(adj_mtx)
# op.to_matrix_op()  # See the operator with matrix form
qp = QuadraticProgram()
qp.from_ising(op)
print(qp)

optimizer = SLSQP(maxiter=10000)
algorithm_globals.random_seed = 1234
seed = 11234
num = 8
backend = Aer.get_backend('statevector_simulator')

ry = TwoLocal(num, 'ry', 'cz', reps=4, entanglement='full')
quantum_instance = QuantumInstance(backend=backend, seed_simulator=seed, seed_transpiler=seed)
vqe = VQE(ry, optimizer=optimizer, quantum_instance=quantum_instance)

vqe_meo = MinimumEigenOptimizer(vqe)

result = vqe_meo.solve(qp)

print(result.samples)

for opt in range(4):
    for it in range(1, 11):
        if opt == 0:
            optimizer = SLSQP(maxiter=1000)  
        elif opt == 1:
            optimizer = COBYLA(maxiter=1000)
        elif opt == 2:
            optimizer = NELDER_MEAD(maxiter=1000)
        elif opt == 3:
            optimizer = POWELL(maxiter=1000)

        algorithm_globals.random_seed = 1234
        seed = 12345
        backend = Aer.get_backend('statevector_simulator')
        qNum = 8
        
        ry = TwoLocal(qNum, 'ry', 'cz', reps=it, entanglement='full')
        quantum_instance = QuantumInstance(backend=backend, seed_simulator=seed, seed_transpiler=seed)
        vqe = VQE(ry, optimizer=optimizer, quantum_instance=quantum_instance)

        vqe_meo = MinimumEigenOptimizer(vqe)

        result = vqe_meo.solve(qp)

        print(result)
        
        # use dictionary form to output result
        output(result, f"data/VQE/prob_all_opt{opt}_layer{it}.npy", f"data/VQE/fval_all_opt{opt}_layer{it}.npy")

classical_exp = 2
probability = np.zeros([4, 10])
cost = np.zeros([4, 10])
pdepth = np.arange(1, 11)
expectation = np.zeros([4, 10])
for opt in range(4):
    for it in range(1, 11):
        df = np.load(f"data/VQE/prob_all_opt{opt}_layer{it}.npy",allow_pickle='TRUE').item()
        df2 = np.load(f"data/VQE/fval_all_opt{opt}_layer{it}.npy",allow_pickle='TRUE').item()
        probability[opt, it-1] = df.get('10011010', 0)
        probability[opt, it-1] += df.get('01100101', 0)
        cost[opt, it-1] = df2.get('10011010', 0)
        for key in df:
            expectation[opt, it-1] += df[key] * df2[key]
        

alpha = (expectation - classical_exp) / np.abs(classical_exp)

plt.figure(figsize=(10,8))
for opt in range(4):
    plt.plot(pdepth, probability[opt])
plt.legend(["SLSQP", "COBYLA", "NELDER_MEAD", "POWELL"])
plt.xlim(0.5, 10.5)
plt.ylim(0, 1)
plt.xticks(np.arange(1, 11), ["1", "2", "3", "4", "5", "6", "7", "8", "9", "10"], fontsize=18)
plt.title("VQE, Probability of Finding Right Answer\n (10011010 or 01100101)", fontsize=18)
plt.xlabel("p-depth", fontsize=18)
plt.ylabel("Probability", fontsize=18)
plt.grid("--")
plt.savefig(f"graph/VQE/prob.png", bbox_inches='tight',pad_inches = 0,dpi=300)
plt.close()

plt.figure(figsize=(10,8))
for opt in range(4):
    plt.plot(pdepth, expectation[opt])
plt.legend(["SLSQP", "COBYLA", "NELDER_MEAD", "POWELL"])
plt.xlim(0.5, 10.5)
plt.ylim(2, 8)
plt.xticks(np.arange(1, 11), ["1", "2", "3", "4", "5", "6", "7", "8", "9", "10"], fontsize=18)
plt.title("VQE, Expectation of Cost Function\n (Right Answer: 2)", fontsize=18)
plt.xlabel("p-depth", fontsize=18)
plt.ylabel("Expectation Value", fontsize=18)
plt.grid("--")
plt.savefig(f"graph/VQE/expect.png", bbox_inches='tight',pad_inches = 0,dpi=300) 
plt.close()

plt.figure(figsize=(15,8))
for opt in range(4):
    plt.plot(pdepth, alpha[opt])
plt.legend(["SLSQP", "COBYLA", "ADAM", "CG", "Gradient descent"])
plt.xlim(0.5, 10.5)
plt.ylim(0, 5)
plt.xticks(np.arange(1, 11), ["1", "2", "3", "4", "5", "6", "7", "8", "9", "10"], fontsize=18)
plt.title(r"VQE, $\frac{expectation - classical\ expectation}{|classical\ expectation|}$", fontsize=18)
plt.xlabel("p-depth", fontsize=18)
plt.ylabel("Error rate", fontsize=18)
plt.grid("--")
plt.savefig(f"graph/VQE/alpha.png", bbox_inches='tight',pad_inches = 0,dpi=300) 
plt.show()
plt.close()

dd = np.load(f"data/QAOA/prob_all_opt0_layer1.npy",allow_pickle='TRUE').item()
new_ans = {}
x = []
for key in dd:
    tmp = 0
    for i in range(8):
        tmp += int(key[i])
    if tmp == 4:
        x.append(key)

for it in range(1,11):
    df = np.load(f"data/VQE/prob_all_opt0_layer{it}.npy",allow_pickle='TRUE').item()
    df2 = np.load(f"data/VQE/prob_all_opt1_layer{it}.npy",allow_pickle='TRUE').item()
    df3 = np.load(f"data/VQE/prob_all_opt2_layer{it}.npy",allow_pickle='TRUE').item()
    df4 = np.load(f"data/VQE/prob_all_opt3_layer{it}.npy",allow_pickle='TRUE').item()

    y = np.zeros(len(x))
    y2 = np.zeros(len(x))
    y3 = np.zeros(len(x))
    y4 = np.zeros(len(x))

    for i in range(len(x)):
        y[i] = float(df.get(x[i], 0.))
        y2[i] = float(df2.get(x[i], 0.))
        y3[i] = float(df3.get(x[i], 0.))
        y4[i] = float(df4.get(x[i], 0.))


    wid = 0.5
    X = np.linspace(0, 150, len(x))
    plt.figure(figsize=(25, 8))
    plt.xlim(-2, 155)
    # plt.ylim(0, 1)
    plt.title(f"VQE, Reps={it}", fontsize=16)
    plt.bar(X, y, width=wid, color='r')
    plt.bar(X+0.5, y2, width=wid, color='b')
    plt.bar(X+1, y3, width=wid, color='k')
    plt.bar(X+1.5, y4, width=wid, color='g')
    plt.legend(labels=["SLSQP", "COBYLA", "NELDER_MEAD", "POWELL"])
    plt.xticks(X, x)
    plt.ylabel("Probabilities", fontsize=16)

    plt.xticks(rotation=90, fontsize=16)
    plt.savefig(f"graph/VQE/few_{it}.png", bbox_inches='tight',pad_inches = 0,dpi=200)
    plt.show()
    plt.close()

lst = [[1,4,6],[0,2],[1,5,7],[4],[0,3,5],[2,4],[0],[2]]
adj_mtx = adjacency_matrix(lst)
op_w = normalized_h(adj_mtx, 1)
# op.to_matrix_op()  # See the operator with matrix form
qp_w = QuadraticProgram()
qp_w.from_ising(op_w)
print(qp_w)

optimizer = SLSQP(maxiter=10000)
algorithm_globals.random_seed = 1234
seed = 11234
num = 8
backend = Aer.get_backend('statevector_simulator')

ry = TwoLocal(num, 'ry', 'cz', reps=4, entanglement='full')
quantum_instance = QuantumInstance(backend=backend, seed_simulator=seed, seed_transpiler=seed)
vqe = VQE(ry, optimizer=optimizer, quantum_instance=quantum_instance)

vqe_meo = MinimumEigenOptimizer(vqe)

result = vqe_meo.solve(qp_w)

print(result.samples)

# VQE data
ratio_arr = np.linspace(0.2,2,5)
for opt in range(2):
    for it in range(1, 11):
        for ratio_idx in range(len(ratio_arr)):
            if opt == 0:
                optimizer = SLSQP(maxiter=1000)  
            elif opt == 1:
                optimizer = COBYLA(maxiter=1000)
            #elif opt == 2:
                #optimizer = NELDER_MEAD(maxiter=1000)
            #elif opt == 3:
             #   optimizer = POWELL(maxiter=1000)

            algorithm_globals.random_seed = 1234
            seed = 12345
            backend = Aer.get_backend('statevector_simulator')
            qNum = 8

            ry = TwoLocal(qNum, 'ry', 'cz', reps=it, entanglement='full')
            quantum_instance = QuantumInstance(backend=backend, seed_simulator=seed, seed_transpiler=seed)
            vqe = VQE(ry, optimizer=optimizer, quantum_instance=quantum_instance)

            vqe_meo = MinimumEigenOptimizer(vqe)
            
            ratio = ratio_arr[ratio_idx]
            op_w = normalized_h(adj_mtx, ratio)
            qp = QuadraticProgram()
            qp.from_ising(op_w)
            result = vqe_meo.solve(qp)

            print(result)

            output(result, f"data_normalized/VQE/prob_all_opt{opt}_layer{it}_ratio{ratio}.npy", f"data_normalized/VQE/fval_all_opt{opt}_layer{it}_ratio{ratio}.npy")

classical_exp = np.load("data/normalized/classical_exp_arr.npy")
classical_exp_arr = np.zeros([50])
for i in range(5):
    classical_exp_arr[10*i:10*(i+1)] = classical_exp[i]
print(classical_exp_arr)

df = np.load(f"data/normalized/VQE/prob_all_opt{0}_layer{1}_ratio{0.2}.npy",allow_pickle='TRUE').item()

classical_exp = np.load("data/normalized/classical_exp_arr.npy")
probability = np.zeros([4, 10])
cost = np.zeros([4, 10])
pdepth = np.arange(1, 11)
expectation = np.zeros([4, 10])

ratMax = [0.2, 0.65, 1.1, 1.55, 2.0]
output_prob = np.zeros([50, 3])   
output_expectation = np.zeros([50, 3])
output_error_rate = np.zeros([50, 3])
count = 0
for opt in range(0, 2):
    for it in range(1, 11):
        for ra in ratMax:
            df = np.load(f"data/normalized/VQE/prob_all_opt{opt}_layer{it}_ratio{ra}.npy",allow_pickle='TRUE').item()
            df2 = np.load(f"data/normalized/VQE/fval_all_opt{opt}_layer{it}_ratio{ra}.npy",allow_pickle='TRUE').item()
            output_prob[count, 0] = output_expectation[count, 0] = output_error_rate[count, 0] = it
            output_prob[count, 1] = output_expectation[count, 1] = output_error_rate[count, 1] = ra
            output_prob[count, 2] = df.get('10011010', 0) + df.get('01100101', 0)
            
            for key in df:
                output_expectation[count, 2] += df[key] * df2[key]
            count += 1

    for i in range(50):
        output_error_rate[i, 2] = (output_expectation[i, 2] - classical_exp_arr[i]) / np.abs(classical_exp_arr[i])
    
    if opt == 0:
        np.savetxt("data/normalized/VQE_output/probability.csv", output_prob, delimiter=",")
        np.savetxt("data/normalized/VQE_output/expectation.csv", output_expectation, delimiter=",")
        np.savetxt("data/normalized/VQE_output/error_rate.csv", output_error_rate, delimiter=",")
        
        np.save("data/normalized/VQE_output_reprocessed/probability.npy", output_prob)
        np.save("data/normalized/VQE_output_reprocessed/expectation.npy", output_expectation)
        np.save("data/normalized/VQE_output_reprocessed/error_rate.npy", output_error_rate)
        
    if opt == 1:
        np.savetxt("data/normalized/VQE_output/probability_2.csv", output_prob, delimiter=",")
        np.savetxt("data/normalized/VQE_output/expectation_2.csv", output_expectation, delimiter=",")
        np.savetxt("data/normalized/VQE_output/error_rate_2.csv", output_error_rate, delimiter=",")
        
        np.save("data/normalized/VQE_output_reprocessed/probability_2.npy", output_prob)
        np.save("data/normalized/VQE_output_reprocessed/expectation_2.npy", output_expectation)
        np.save("data/normalized/VQE_output_reprocessed/error_rate_2.npy", output_error_rate)

lst = [[1,4,6],[0,2],[1,5,7],[4],[0,3,5],[2,4],[0],[2]]
adj_mtx = adjacency_matrix(lst)
op = unnormalized_h(adj_mtx)
# op.to_matrix_op()  # See the operator with matrix form
qp = QuadraticProgram()
qp.from_ising(op)
print(qp)

optimizer = SLSQP(maxiter=10000)
algorithm_globals.random_seed = 1234
seed = 11234
backend = Aer.get_backend('statevector_simulator')

quantum_instance = QuantumInstance(backend=backend, seed_simulator=seed, seed_transpiler=seed)
qaoa = QAOA(optimizer=optimizer, reps=3, quantum_instance=quantum_instance)
meo = MinimumEigenOptimizer(qaoa)

qaoa_meo = MinimumEigenOptimizer(qaoa)

result = qaoa_meo.solve(qp)

print(result.samples)

for opt in range(4):
    for it in range(1, 11):
        if opt == 0:
            optimizer = SLSQP(maxiter=1000)  
        elif opt == 1:
            optimizer = COBYLA(maxiter=1000)
        elif opt == 2:
            optimizer = NELDER_MEAD(maxiter=1000)
        elif opt == 3:
            optimizer = POWELL(maxiter=1000)

        algorithm_globals.random_seed = 1234
        seed = 12345
        backend = Aer.get_backend('statevector_simulator')
        qNum = 8
        
        qaoa = QAOA(optimizer=optimizer, reps=3, quantum_instance=quantum_instance)
        meo = MinimumEigenOptimizer(qaoa)
        qaoa_meo = MinimumEigenOptimizer(qaoa)

        result = qaoa_meo.solve(qp)

        print(result)

        output(result, f"data/QAOA/prob_all_opt{opt}_layer{it}.npy", f"data/QAOA/fval_all_opt{opt}_layer{it}.npy")

classical_exp = 2
probability = np.zeros([4, 10])
cost = np.zeros([4, 10])
pdepth = np.arange(1, 11)
expectation = np.zeros([4, 10])
for opt in range(4):
    for it in range(1, 11):
        df = np.load(f"data/QAOA/prob_all_opt{opt}_layer{it}.npy",allow_pickle='TRUE').item()
        df2 = np.load(f"data/QAOA/fval_all_opt{opt}_layer{it}.npy",allow_pickle='TRUE').item()
        probability[opt, it-1] = df.get('10011010', 0)
        probability[opt, it-1] += df.get('01100101', 0)
        cost[opt, it-1] = df2.get('10011010', 0)
        for key in df:
            expectation[opt, it-1] += df[key] * df2[key]
        

alpha = (expectation - classical_exp) / np.abs(classical_exp)

plt.figure(figsize=(10,8))
for opt in range(4):
    plt.plot(pdepth, probability[opt])
plt.legend(["SLSQP", "COBYLA", "NELDER_MEAD", "POWELL"])
plt.xlim(0.5, 10.5)
plt.ylim(0, 1)
plt.xticks(np.arange(1, 11), ["1", "2", "3", "4", "5", "6", "7", "8", "9", "10"], fontsize=18)
plt.title("QAOA, Probability of Finding Right Answer\n (10011010 or 01100101)", fontsize=18)
plt.xlabel("p-depth", fontsize=18)
plt.ylabel("Probability", fontsize=18)
plt.grid("--")
plt.savefig(f"graph/QAOA/prob.png", bbox_inches='tight',pad_inches = 0,dpi=300)
plt.close()

plt.figure(figsize=(10,8))
for opt in range(4):
    plt.plot(pdepth, expectation[opt])
plt.legend(["SLSQP", "COBYLA", "NELDER_MEAD", "POWELL"])
plt.xlim(0.5, 10.5)
# plt.ylim(2, 8)
plt.xticks(np.arange(1, 11), ["1", "2", "3", "4", "5", "6", "7", "8", "9", "10"], fontsize=18)
plt.title("QAOA, Expectation of Cost Function\n (Right Answer: 2)", fontsize=18)
plt.xlabel("p-depth", fontsize=18)
plt.ylabel("Expectation Value", fontsize=18)
plt.grid("--")
plt.savefig(f"graph/QAOA/expect.png", bbox_inches='tight',pad_inches = 0,dpi=300) 
plt.close()

plt.figure(figsize=(15,8))
for opt in range(4):
    plt.plot(pdepth, alpha[opt])
plt.legend(["SLSQP", "COBYLA", "NELDER_MEAD", "POWELL"])
plt.xlim(0.5, 10.5)
plt.ylim(0, 5)
plt.xticks(np.arange(1, 11), ["1", "2", "3", "4", "5", "6", "7", "8", "9", "10"], fontsize=18)
plt.title(r"QAOA, $\frac{expectation - classical\ expectation}{|classical\ expectation|}$", fontsize=18)
plt.xlabel("p-depth", fontsize=18)
plt.ylabel("Error rate", fontsize=18)
plt.grid("--")
plt.savefig(f"graph/QAOA/alpha.png", bbox_inches='tight',pad_inches = 0,dpi=300) 
plt.show()
plt.close()

lst = [[1,4,6],[0,2],[1,5,7],[4],[0,3,5],[2,4],[0],[2]]
adj_mtx = adjacency_matrix(lst)
op_w = normalized_h(adj_mtx, 1)
# op.to_matrix_op()  # See the operator with matrix form
qp_w = QuadraticProgram()
qp_w.from_ising(op_w)
print(qp_w)

optimizer = SLSQP(maxiter=10000)
algorithm_globals.random_seed = 1234
seed = 11234
backend = Aer.get_backend('statevector_simulator')

quantum_instance = QuantumInstance(backend=backend, seed_simulator=seed, seed_transpiler=seed)
qaoa = QAOA(optimizer=optimizer, reps=3, quantum_instance=quantum_instance)
meo = MinimumEigenOptimizer(qaoa)

qaoa_meo = MinimumEigenOptimizer(qaoa)

result = qaoa_meo.solve(qp_w)
print(result.samples)

x = []
prob = np.zeros([4, 10])
func = np.zeros([4, 10])
ratio_arr = np.linspace(0.2,2,5)

for opt in range(2):
    for it in range(1, 11):
        for ratio_idx in range(len(ratio_arr)):
            if opt == 0:
                optimizer = SLSQP(maxiter=1000)  
            elif opt == 1:
                optimizer = COBYLA(maxiter=1000)
            # if opt == 2:
            #    optimizer = NELDER_MEAD(maxiter=1000)
            # elif opt == 3:
            #    optimizer = POWELL(maxiter=1000)

            algorithm_globals.random_seed = 1234
            seed = 12345
            backend = Aer.get_backend('statevector_simulator')
            qNum = 8
            
            quantum_instance = QuantumInstance(backend=backend, seed_simulator=seed, seed_transpiler=seed)
            qaoa = QAOA(optimizer=optimizer, reps=it, quantum_instance=quantum_instance)
            meo = MinimumEigenOptimizer(qaoa)
            qaoa_meo = MinimumEigenOptimizer(qaoa)

            #prepare qp
            ratio = ratio_arr[ratio_idx]
            op_w = normalized_h(adj_mtx, ratio)
            qp = QuadraticProgram()
            qp.from_ising(op_w)
            result = qaoa_meo.solve(qp)

            print(result)
            
            output(result, f"data_normalized/QAOA/prob_all_opt{opt}_layer{it}_ratio{ratio}.npy", f"data_normalized/QAOA/fval_all_opt{opt}_layer{it}_ratio{ratio}.npy")

df = np.load(f"data/normalized/QAOA/prob_all_opt{0}_layer{1}_ratio{0.2}.npy",allow_pickle='TRUE').item()

classical_exp = np.load("data/normalized/classical_exp_arr.npy")
probability = np.zeros([4, 10])
cost = np.zeros([4, 10])
pdepth = np.arange(1, 11)
expectation = np.zeros([4, 10])

ratMax = [0.2, 0.65, 1.1, 1.55, 2.0]
output_prob = np.zeros([50, 3])   
output_expectation = np.zeros([50, 3])
output_error_rate = np.zeros([50, 3])
count = 0
for opt in range(0, 1):
    for it in range(1, 11):
        for ra in ratMax:
            df = np.load(f"data/normalized/QAOA/prob_all_opt{opt}_layer{it}_ratio{ra}.npy",allow_pickle='TRUE').item()
            df2 = np.load(f"data/normalized/QAOA/fval_all_opt{opt}_layer{it}_ratio{ra}.npy",allow_pickle='TRUE').item()
            output_prob[count, 0] = output_expectation[count, 0] = output_error_rate[count, 0] = it
            output_prob[count, 1] = output_expectation[count, 1] = output_error_rate[count, 1] = ra
            output_prob[count, 2] = df.get('10011010', 0) + df.get('01100101', 0)
            
            for key in df:
                output_expectation[count, 2] += df[key] * df2[key]
            count += 1

    for i in range(50):
        output_error_rate[i, 2] = (output_expectation[i, 2] - classical_exp_arr[i]) / np.abs(classical_exp_arr[i])

    if opt == 0:
        np.savetxt("data/normalized/QAOA_output/probability.csv", output_prob, delimiter=",")
        np.savetxt("data/normalized/QAOA_output/expectation.csv", output_expectation, delimiter=",")
        np.savetxt("data/normalized/QAOA_output/error_rate.csv", output_error_rate, delimiter=",")
        
        np.save("data/normalized/QAOA_output_reprocessed/probability.npy", output_prob)
        np.save("data/normalized/QAOA_output_reprocessed/expectation.npy", output_expectation)
        np.save("data/normalized/QAOA_output_reprocessed/error_rate.npy", output_error_rate)
    
    
    if opt == 1:
        np.savetxt("data/normalized/QAOA_output/probability_2.csv", output_prob, delimiter=",")
        np.savetxt("data/normalized/QAOA_output/expectation_2.csv", output_expectation, delimiter=",")
        np.savetxt("data/normalized/QAOA_output/error_rate_2.csv", output_error_rate, delimiter=",")
        
        np.save("data/normalized/QAOA_output_reprocessed/probability_2.npy", output_prob)
        np.save("data/normalized/QAOA_output_reprocessed/expectation_2.npy", output_expectation)
        np.save("data/normalized/QAOA_output_reprocessed/error_rate_2.npy", output_error_rate)

iter_arr = np.linspace(1,8,8)
ratio_arr = np.linspace(0.2,2,5)
op = operator_from_adjacency_matrix(adj_mtx)
for iter_idx in range(len(iter_arr)):
        
    iter_num = iter_arr[iter_idx]
    qp = QuadraticProgram()
    qp.from_ising(op)
    backend = Aer.get_backend('statevector_simulator')

    grover_optimizer = GroverOptimizer(8, num_iterations=iter_num, quantum_instance=backend)
    result = grover_optimizer.solve(qp)
    print(result.prettyprint())
    
    output(result, f"data/GROVER/prob_iter{iter_num}.npy", f"data/GROVER/fval_iter{iter_num}.npy")

plt.figure(figsize=(15,8))
legend_arr = []
for i in range(len(ratio_arr)):
    legend_arr.append("ratio = "+str(ratio_arr[i]))
    plt.plot(iter_arr, alpha[i])
legend_arr.append("Unnormalized")
plt.plot(iter_arr,unnormalized_alpha)

xticks_arr = []
for i in range(len(iter_arr)):
    xticks_arr.append(str(iter_arr[i]))
    
plt.legend(legend_arr)
plt.xlim(0.5, 8.5)
plt.xticks(iter_arr, xticks_arr, fontsize=18)
plt.title(r"GROVER, $\frac{expectation - classical\ expectation}{|classical\ expectation|}$", fontsize=18)
plt.xlabel("Number of iteration", fontsize=18)
plt.ylabel("Error rate", fontsize=18)
plt.grid("--")
plt.savefig(f"graph/GROVER/NORMALIZED/error_rate.png", bbox_inches='tight',pad_inches = 0,dpi=300)
plt.show()

normalized_error_rate = np.load("data/normalized/VQE_output_reprocessed/error_rate.npy")
# df = pd.read_csv('data/normalized/VQE_output/probability.csv', sep=',',header=None)
x = np.arange(1, 11)
y = np.zeros([5, 10])
ratio_arr = np.linspace(0.2,2,5)
# print(normalized_error_rate)

normalized_error_rate = normalized_error_rate[:, 2:]
for i in range(5):
    y[i] = normalized_error_rate[i::5].T


    
    
classical_exp = 2
probability = np.zeros([4, 10])
cost = np.zeros([4, 10])
pdepth = np.arange(1, 11)
expectation = np.zeros([4, 10])
for opt in range(4):
    for it in range(1, 11):
        df = np.load(f"data/VQE/prob_all_opt{opt}_layer{it}.npy",allow_pickle='TRUE').item()
        df2 = np.load(f"data/VQE/fval_all_opt{opt}_layer{it}.npy",allow_pickle='TRUE').item()
        probability[opt, it-1] = df.get('10011010', 0)
        probability[opt, it-1] += df.get('01100101', 0)
        cost[opt, it-1] = df2.get('10011010', 0)
        for key in df:
            expectation[opt, it-1] += df[key] * df2[key]
        

alpha = (expectation - classical_exp) / np.abs(classical_exp)

plt.figure(figsize=(15,8))
for i in range(5):
    plt.plot(x, y[i])
plt.plot(x, alpha[0])
plt.legend(["ratio = 0.2", "ratio = 0.65", "ratio = 1.1", "ratio = 1.55", "ratio = 2", "Unnormalized"])
plt.xlim(0.5, 10.5)
# plt.ylim(0, 5)
plt.xticks(np.arange(1, 11), ["1", "2", "3", "4", "5", "6", "7", "8", "9", "10"], fontsize=18)
plt.title(r"VQE with SLSQP, $\frac{expectation - classical\ expectation}{|classical\ expectation|}$", fontsize=18)
plt.xlabel("p-depth", fontsize=18)
plt.ylabel("Error rate", fontsize=18)
plt.grid("--")
plt.savefig(f"graph/VQE_compare.png", bbox_inches='tight',pad_inches = 0,dpi=300) 
plt.show()
plt.close()


normalized_error_rate = np.load("data/normalized/QAOA_output_reprocessed/error_rate.npy")
# df = pd.read_csv('data/normalized/VQE_output/probability.csv', sep=',',header=None)
x = np.arange(1, 11)
y = np.zeros([5, 10])
ratio_arr = np.linspace(0.2,2,5)
# print(normalized_error_rate)

normalized_error_rate = normalized_error_rate[:, 2:]
for i in range(5):
    y[i] = normalized_error_rate[i::5].T

    
classical_exp = 2
probability = np.zeros([4, 10])
cost = np.zeros([4, 10])
pdepth = np.arange(1, 11)
expectation = np.zeros([4, 10])
for opt in range(4):
    for it in range(1, 11):
        df = np.load(f"data/QAOA/prob_all_opt{opt}_layer{it}.npy",allow_pickle='TRUE').item()
        df2 = np.load(f"data/QAOA/fval_all_opt{opt}_layer{it}.npy",allow_pickle='TRUE').item()
        probability[opt, it-1] = df.get('10011010', 0)
        probability[opt, it-1] += df.get('01100101', 0)
        cost[opt, it-1] = df2.get('10011010', 0)
        for key in df:
            expectation[opt, it-1] += df[key] * df2[key]
        

alpha = (expectation - classical_exp) / np.abs(classical_exp)

plt.figure(figsize=(15,8))
for i in range(5):
    plt.plot(x, y[i])
plt.plot(x, alpha[0])
plt.legend(["ratio = 0.2", "ratio = 0.65", "ratio = 1.1", "ratio = 1.55", "ratio = 2", "Unnormalized"])
plt.xlim(0.5, 10.5)
# plt.ylim(0, 5)
plt.xticks(np.arange(1, 11), ["1", "2", "3", "4", "5", "6", "7", "8", "9", "10"], fontsize=18)
plt.title(r"QAOA with SLSQP, $\frac{expectation - classical\ expectation}{|classical\ expectation|}$", fontsize=18)
plt.xlabel("p-depth", fontsize=18)
plt.ylabel("Error rate", fontsize=18)
plt.grid("--")
plt.savefig(f"graph/QAOA_compare.png", bbox_inches='tight',pad_inches = 0,dpi=300) 
plt.show()
plt.close()


n = 8
penalty = 2 * n

linear2penalty = LinearEqualityToPenalty(penalty=penalty)
qp = linear2penalty.convert(qp)
_, offset = qp.to_ising()

# set classical optimizer
maxiter = 1000
optimizer = COBYLA(maxiter=maxiter)
# set variational ansatz
ansatz = RealAmplitudes(n, reps=1)
m = ansatz.num_parameters

# set backend
backend_name = "qasm_simulator"  # use this for QASM simulator
# backend_name = 'aer_simulator_statevector'  # use this for statevector simlator
backend = Aer.get_backend(backend_name)

# run variational optimization for different values of alpha
alphas = [1.0, 0.50, 0.25]  # confidence levels to be evaluated
# dictionaries to store optimization progress and results
objectives = {alpha: [] for alpha in alphas}  # set of tested objective functions w.r.t. alpha
results = {}  # results of minimum eigensolver w.r.t alpha

# callback to store intermediate results
def callback(i, params, obj, stddev, alpha):
    # we translate the objective from the internal Ising representation
    # to the original optimization problem
    objectives[alpha] += [-(obj + offset)]

for alpha in alphas:

    # initialize CVaR_alpha objective
    cvar_exp = CVaRExpectation(alpha, PauliExpectation())
    cvar_exp.compute_variance = lambda x: [0]  # to be fixed in PR #1373

    # initialize VQE using CVaR
    vqe = VQE(
        expectation=cvar_exp,
        optimizer=optimizer,
        ansatz=ansatz,
        quantum_instance=backend,
        callback=lambda i, params, obj, stddev: callback(i, params, obj, stddev, alpha),
    )

    # initialize optimization algorithm based on CVaR-VQE
    opt_alg = MinimumEigenOptimizer(vqe)

    # solve problem
    results[alpha] = opt_alg.solve(qp)

    # print results
    print("alpha = {}:".format(alpha))
    print(results[alpha].prettyprint())
    print()