import numpy as np

from qiskit import Aer, IBMQ
from qiskit.aqua import aqua_globals, QuantumInstance
from qiskit.aqua.algorithms import QAOA
from qiskit.aqua.components.optimizers import *
from qiskit.quantum_info import Pauli
from qiskit.aqua.operators import WeightedPauliOperator
from qiskit.providers.aer.noise import NoiseModel

provider = IBMQ.load_account()

def get_operator(J,h,n):    
    pauli_list = []

    for (i,j) in J: # For each coefficient in J (couplings) we add a term J[i,j]Z_iZj
        x_p = np.zeros(n, dtype=np.bool)
        z_p = np.zeros(n, dtype=np.bool)
        z_p[n-1-i] = True 
        z_p[n-1-j] = True
        pauli_list.append([J[(i,j)],Pauli(z_p, x_p)])
     
    for i in h: # For each coefficient in h we add a term h[i]Z_i
        x_p = np.zeros(n, dtype=np.bool)
        z_p = np.zeros(n, dtype=np.bool)
        z_p[n-1-i] = True
        pauli_list.append([h[i],Pauli(z_p, x_p)])
    
    return WeightedPauliOperator(paulis=pauli_list)

# Edges of the graph

J1 = {(0,1):1, (1,2):1, (2,3):1, (3,4):1, (4,0):1}
h1 = {}
n = 5

# Hamiltonian

q_op =get_operator(J1,h1,n) 
print(q_op)
q_op.print_details()

rep = 10
backend = Aer.get_backend('statevector_simulator')
quantum_instance = QuantumInstance(backend)

p = 1
val = 0
for i in range(rep):
    print("----- ITERATION ",i, " ------")
    optimizer = COBYLA()
    qaoa = QAOA(q_op, optimizer, p=p)
    result = qaoa.run(quantum_instance)
    print("Optimal value", result['optimal_value'])
    val+=result['optimal_value']
print("----- AVERAGE -----")
print("Average value",val/rep)

p = 2
val = 0
for i in range(rep):
    print("----- ITERATION ",i, " ------")
    optimizer = COBYLA()
    qaoa = QAOA(q_op, optimizer, p=p)
    result = qaoa.run(quantum_instance)
    print("Optimal value", result['optimal_value'])
    val+=result['optimal_value']
print("----- AVERAGE -----")
print("Average value",val/rep)

rep = 10
backendIBM = provider.get_backend('ibmq_ourense')
noise_model = NoiseModel.from_backend(backendIBM)
coupling_map = backendIBM.configuration().coupling_map
basis_gates = noise_model.basis_gates
backend = Aer.get_backend("qasm_simulator")


shots = 8192
optimization_level = 3
p = 1
quantum_instance = QuantumInstance(backend, shots = shots, 
                                    optimization_level = optimization_level,
                                    noise_model = noise_model,
                                    basis_gates = basis_gates,
                                    coupling_map = coupling_map)

p = 1
val = 0
for i in range(rep):
    print("----- ITERATION ",i, " ------")
    optimizer = COBYLA()
    qaoa = QAOA(q_op, optimizer, p=p)
    result = qaoa.run(quantum_instance)
    print("Optimal value", result['optimal_value'])
    val+=result['optimal_value']
print("----- AVERAGE -----")
print("Average value",val/rep)