%config IPCompleter.greedy=True
%matplotlib inline

from qplex_core import *

precision = 6
coefs_param = [2, -3]
beta = 7
a = wrapper_optimiza_f(precision, coefs_param, beta)
print(a)

#tú mandas la función 2*X-3*Y,
#La precisión para representar enteros es de 6 qbits por número (usamos los X e Y y le añadimos un Z)
#como restricción se manda X*y<=6
#El resultado es {0: 0, 1: 7}, que quiere decir, que X vale 0 e Y vale 7




# useful additional packages 
import matplotlib.pyplot as plt
import matplotlib.axes as axes
import numpy as np
import networkx as nx

from qiskit import BasicAer
from qiskit.tools.visualization import plot_histogram
from qiskit.aqua import Operator, run_algorithm
from qiskit.aqua.input import EnergyInput
from qiskit.aqua.translators.ising import max_cut, tsp
from qiskit.aqua.algorithms import VQE, ExactEigensolver
from qiskit.aqua.components.optimizers import SPSA
from qiskit.aqua.components.variational_forms import RY
from qiskit.aqua import QuantumInstance

# setup aqua logging
import logging
from qiskit.aqua import set_qiskit_aqua_logging
# set_qiskit_aqua_logging(logging.DEBUG)  # choose INFO, DEBUG to see the log

from docplex.mp.model import Model
from qiskit.aqua.translators.ising import docplex

from qplex_core import *

def mainJGOS():


    precision = 6
    
    coefs_param = [2, -3 ]
    coefs = coefs_param.copy()
    coefs.append(0)
    coefs_restr = (1, 1, 1)
    beta = 7
    lista_vars = ['a', 'b', 'c', 'd', 'e','f','g','h']
    mdl = Model()
    variables = []
    for num_var in range(len(coefs)):
        tmp ={i: mdl.binary_var(name=(lista_vars[num_var]+'_{0}').format(i)) for i in range(precision)}
        variables.append(tmp)
    #x = {i: mdl.binary_var(name='x_{0}'.format(i)) for i in range(precision)}
    #y = {i: mdl.binary_var(name='y_{0}'.format(i)) for i in range(precision)}
    #z = {i: mdl.binary_var(name='z_{0}'.format(i)) for i in range(precision)}
    print(variables)
    # Object function
    #my_func = mdl.sum(coefs[0]*(2**i)*x[i]+coefs[1]*(2**i)*y[i]+(2**i)*coefs[2]*z[i] for i in range(precision))

    #my_func = mdl.sum(coefs[j]*(2**i)*vars[j][i] for j in range(len(coefs)) for i in range(precision))
    
    
    my_func = mdl.sum(coefs[j]*(2**i)*variables[j][i] for j in range(len(coefs)) for i in range(precision))
    
    
    #(x[i] for i in range(precision)), (y[i] for i in range(precision)), (z[i] for i in range(precision)
    #tmp = {0:{'var':x,'coef':coefs_restr[0]}, 
    #       1:{'var':y,'coef':coefs_restr[1]}, 
    #       2:{'var':z,'coef':coefs_restr[2]}}                          
    
    mdl.maximize(my_func)


    inverted_beta = dameInversoBinario(beta, precision, len(coefs))
    
    #mdl.add_constraint(mdl.sum( tmp[v]['var'][i]*(2**i)*tmp[v]['coef'] for v in range(len(coefs)) for i in range(precision)) == beta)
    #mdl.add_constraint(mdl.sum( x[0] + x[1]*(2) + x[2]*(4) + y[0] + y[1]*(2) + y[2]*(4) + z[0] + z[1]*(2) + z[2]*(4) ) == inverted_beta)
    #mdl.add_constraint(mdl.sum( variables[0][0] + variables[0][1]*(2) + variables[0][2]*(4) + variables[1][0] + variables[1][1]
    #*(2) + variables[1][2]*(4) + variables[2][0] + variables[2][1]*(2) + variables[2][2]*(4) ) == inverted_beta)
    
    mdl.add_constraint(mdl.sum(variables[v][i]*2**i  for v in range(len(coefs)) for i in range(precision)) == inverted_beta)
    
    #mdl.add_constraint(mdl.sum( -x[0] - x[1]*(2) - x[2]*(4) - y[0] - y[1]*(2) - y[2]*(4) - z[0] - z[1]*(2) - z[2]*(4) ) == 6)

    #mdl.add_constraint(mdl.sum( -1*tmp[v]['var'][i]*(2**i)*tmp[v]['coef'] for v in range(len(coefs)) for i in range(precision)) == -beta)

    qubitOp_docplex, offset_docplex = docplex.get_qubitops(mdl)

    print(qubitOp_docplex)

    #algo_input = EnergyInput(qubitOp_docplex)
    #print(algo_input.)


    #ee = VQE(qubitOp_docplex)
    #ee.run()
    ee = ExactEigensolver(qubitOp_docplex, k=1)
    result_ee = ee.run()
    x_ee = max_cut.sample_most_likely(result_ee['eigvecs'][0])
    print('solution:', max_cut.get_graph_solution(x_ee))
    solucion_ee =  max_cut.get_graph_solution(x_ee)
    return (solucion_ee, None)
    """
    algorithm_cfg = {
        'name': 'ExactEigensolver',
    }
    
    params = {
        'problem': {'name': 'ising'},
        'algorithm': algorithm_cfg
    }
    result = run_algorithm(params,algo_input)
    """
    #x = max_cut.sample_most_likely(result['eigvecs'][0])
    #print('energy:', result['energy'])
    #print('max-cut objective:', result['energy'] + offset_docplex)
    #print('solution:', max_cut.get_graph_solution(x))
    #print('solution objective:', max_cut.max_cut_value(x, w))

    seed = 10598
    
    #change optimizer(spsa), change ry (riyc)
    spsa = SPSA(max_trials=300)
    ry = RY(qubitOp_docplex.num_qubits, depth=6, entanglement='linear')
    vqe = VQE(qubitOp_docplex, ry, spsa, 'matrix')

    backend = BasicAer.get_backend('statevector_simulator')
    quantum_instance = QuantumInstance(backend, seed=seed, seed_transpiler=seed)

    result = vqe.run(quantum_instance)
    x = max_cut.sample_most_likely(result['eigvecs'][0])
    print('solution:', max_cut.get_graph_solution(x))
    return (solucion_ee, max_cut.get_graph_solution(x))

resultado = mainJGOS()
r_ee = resultado[0]
vars = np.array_split(r_ee, 3)
tam = len(vars)
for i in range(tam):
    curr = vars[i]
    curr_i = np.flip(curr)
    curr_i = [ str(int(k)) for k in curr_i ]
    en_binario = "".join(list(curr_i))
    curr_var = int(en_binario, 2)
    if i != tam-1:
        print("Var %d = %d" %(i,curr_var))