from functions import BinPacking, BinPackingNewApproach, new_eq_optimal,  get_figure, interpret, eval_constrains
from functions import mapping_cost, cost_func, qaoa_circuit, check_best_sol
from qiskit.algorithms import QAOA, NumPyMinimumEigensolver
from qiskit_optimization.algorithms import CplexOptimizer, MinimumEigenOptimizer
from qiskit.algorithms.optimizers import COBYLA
import numpy as np
from qiskit_optimization.problems.constraint import ConstraintSense
import matplotlib.pyplot as plt
from scipy.optimize import minimize

from qiskit import Aer
backend = Aer.get_backend("qasm_simulator")

#np.random.seed(1)
num_items = 4 # number of items
num_bins = num_items # maximum number of bins
max_weight = 15 # max weight of a bin
cases = 1
solutions_new = {}
optimal_new = []
ratio_new = []; new =[]
result_classical = []; result_qaoa_new = []
check_const = []
weights=[]
qaoa = MinimumEigenOptimizer(QAOA(optimizer=COBYLA(maxiter=100), reps=1, quantum_instance=backend))
for i in range(cases): # Testing 5 different randomly selected configurations of the problem
    print(f"----------- Case {i+1} -------------")
    weights.append(np.random.randint(1, max_weight, num_items)) # Randomly picking the item weight
    qubo_new = BinPackingNewApproach(num_items, num_bins, weights[-1], max_weight, alpha=1, simplification=True)
    qubo_classical, qp = BinPacking(num_items, num_bins, weights[-1], max_weight, simplification=True)
    result_classical.append(CplexOptimizer().solve(qubo_classical))
    optimal_new.append(new_eq_optimal(qubo_new, qubo_classical))
    result_qaoa_new.append(qaoa.solve(qubo_new))
    solutions_new = result_qaoa_new[-1].fval
    new.append(solutions_new)
    check_const.append(eval_constrains(qp, result_qaoa_new[-1].x, max_weight))
    print(check_const[-1])


print(f" The percetage of positive cases is {100*np.sum(check_const)/cases} %")

for i in range(len(result_classical)):
    fig = get_figure(interpret(result_classical[i].x, weights[i], max_weight, num_items, num_bins, simplify=True), weights[i], max_weight, title=f"-------Case {i+1}-----")
    fig = get_figure(interpret(result_qaoa_new[i].x, weights[i], max_weight, num_items, num_bins, simplify=True), weights[i], max_weight, title="New")

n = 50
alpha = np.linspace(0, 2*np.pi, n)
beta = np.linspace(0, np.pi, n)
map_cost_10 = mapping_cost(alpha, beta, qubo_new, n = 10)
map_cost_20 = mapping_cost(alpha, beta, qubo_new, n = 20)

fig, ax = plt.subplots(1,2, figsize=(10,10))
ax[0].imshow(np.log(map_cost_10))
ax[1].imshow(np.log(map_cost_20))
for i in range(2):
#     ax[i].set_xticks([0,9,19], ["0", r"$\pi/2$", r"$\pi$"])
#     ax[i].set_yticks([0,9,19], ["0", r"$\pi$", r"2$\pi$"])
    ax[i].set_xlabel(r"$\beta$", fontsize=22)
    ax[i].set_ylabel(r"$\alpha$", fontsize=22)


p = 4
x0 = np.random.rand(p*2)
circuit = qaoa_circuit(qubo_new, p=p)
objective = qubo_new.objective
num_evals = 20
fmin = minimize(cost_func, x0, args = (circuit, objective, num_evals), method="COBYLA")

best_sol = check_best_sol(fmin.x, circuit, qp, max_weight, n=10)

print(f" minimum my cost_fun {best_sol}, minimum with cplex {result_qaoa_new[-1].x}")

fig = get_figure(interpret(np.array(best_sol), weights[-1], max_weight, num_items, num_bins, simplify=True), weights[-1], max_weight, title=f"-------Case {i+1}-----")


best_sol