import numpy as np
from qiskit_optimization import QuadraticProgram
from qiskit_optimization.converters import QuadraticProgramToQubo
import dimod

qp = QuadraticProgram("Example")
n = 3
for j in range(n):
            qp.binary_var('x' + "_" + str(j))
linear_objective = np.array([-1,0,-1])
Q = np.array([[-2,0,1],
              [0,2,-2],
              [1,-2,0]]) / 2
quadratic_objective = {}
for i in range(n):
    for j in range(i, n):
        var_i = 'x' + "_" + str(i)
        if i == j:
            quadratic_objective[var_i, var_i] = Q[i, i]
        else:
            var_j = 'x' + "_" + str(j)
            quadratic_objective[var_i, var_j] = Q[i, j] + Q[j, i]
qp.linear_constraint(linear=np.array([0,1,-1]), sense='==', rhs=0,
                        name='linear_cstr')
qp.minimize(linear=linear_objective, quadratic=quadratic_objective)
print(qp) 

conv = QuadraticProgramToQubo()
qubo = conv.convert(qp)
print(f"Auto calculated {conv.penalty=}")
print(qubo)

bqm_binary = dimod.as_bqm(qubo.objective.linear.to_array(), qubo.objective.quadratic.to_array(), dimod.BINARY)
print(bqm_binary)
bqm_ising = bqm_binary.to_ising()
print("Ising Model form coefficients:")
print(bqm_ising)