import numpy as np
w = np.array([[0., 1., 1., 1.], 
     [1., 0., 1., 0.],
     [1., 1., 0., 1.],
     [1., 0., 1., 0.]])

from qiskit_optimization.applications import Maxcut
max_cut = Maxcut(w)
qp = max_cut.to_quadratic_program()
print(qp.export_as_lp_string())

qubitOp, offset = qp.to_ising()
num_qubits = qubitOp.num_qubits
print("Offset:", offset)
print("Ising Hamiltonian:")
print(str(qubitOp))

from qiskit_optimization.algorithms import MinimumEigenOptimizer
from qiskit.algorithms import NumPyMinimumEigensolver
exact = MinimumEigenOptimizer(NumPyMinimumEigensolver())
result = exact.solve(qp)
print(result)

from qiskit.opflow import PauliSumOp
import numpy as np

ops = []
h = qubitOp.primitive
h.coeffs = np.real(qubitOp.coeffs)
ops.append(PauliSumOp(h))
print(ops)

from qiskit_ibm_runtime import QiskitRuntimeService
service = QiskitRuntimeService()
options = {"backend" : "ibmq_qasm_simulator"}

from qiskit.circuit.library import TwoLocal
circuit = TwoLocal(qubitOp.num_qubits, "ry", "cz", reps=5, entanglement="linear")
circuit.decompose().draw(fold=300)

from qiskit_ibm_runtime import Estimator
from qiskit.algorithms.optimizers import SPSA

history = {"nfevs": [], "points": [], "fvals": [], "updates": [], "accepted": []}

with Estimator(
    circuits=circuit,
    observables=ops,
    options=options,
    service=service
) as estimator:
    
    def evaluate_expectation(x):
        x = list(x)
        results = estimator(circuits=[0], observables=[0], parameter_values=[x]).values[0]
        return np.real(results)

    def callback(nfev, point, fval, update, accepted):
        print('expectation value: {}'.format(fval))
        history["nfevs"].append(nfev)
        history["points"].append(point)
        history["fvals"].append(fval)
        history["updates"].append(update)
        history["accepted"].append(accepted)    

    np.random.seed(1)
    initial_point = np.random.uniform(-np.pi, np.pi, len(circuit.parameters))
    optimizer = SPSA(50, callback=callback)
    result = optimizer.minimize(fun=evaluate_expectation, x0=initial_point)


print('optimized circuit parameters:\n{}'.format(result.x))

from qiskit_ibm_runtime import Sampler

circut = circuit.measure_all()
with Sampler(
    circuits=circuit,
    options=options,
    service=service
) as sampler:
    sampler_result = sampler(circuits=[0], parameter_values=[result.x], shots=8000)
    print(result)

bitstrings = sampler_result.quasi_dists[0] # get the dictionary of bitstring counts

import matplotlib.pyplot as plt

fig, axs = plt.subplots(1, 2, figsize=(14, 5))
axs[0].plot(history['fvals'])
axs[1].bar(list(bitstrings.keys()), list(bitstrings.values()))
axs[0].set_xlabel("Energy Evaluation Number")
axs[0].set_ylabel("Energy")
axs[1].set_xlabel("Objective Function Value")
axs[1].set_ylabel("Probability")

samples = []
for x, probability in bitstrings.items():

    x = list(map(int, x)) # convert bitstring to list
    status = qp.is_feasible(x)
    fval = qp.objective.evaluate(x)

    samples.append((x, fval, probability, status))

sorted_samples = sorted(
    samples,
    key=lambda v: (1-v[3], qp.objective.sense.value * v[1]) # feasible: 0; infeasible: 1
)

print('\nsorted')
for sample in sorted_samples:
    print(sample)

best = sorted_samples[0]
print('optimal value: {}'.format(best[0]))
print('optimal function value: {}'.format(best[1]))
print('probability of optimal value: {}'.format(best[2]))

from qiskit_optimization import QuadraticProgram

qp = QuadraticProgram()
qp.binary_var("x")
qp.binary_var("y")
qp.binary_var("z")
qp.minimize(linear=[1, -2, 3], quadratic={("x", "y"): 1, ("x", "z"): -1, ("y", "z"): 2})
print(qp.export_as_lp_string())

qubitOp, offset = qp.to_ising()
num_qubits = qubitOp.num_qubits
print("offset: {}".format(offset))
print("operator:\n{}".format(qubitOp))

from qiskit_optimization.algorithms import MinimumEigenOptimizer
from qiskit.algorithms import NumPyMinimumEigensolver
exact = MinimumEigenOptimizer(NumPyMinimumEigensolver())
result = exact.solve(qp)
print(result)

from qiskit.opflow import PauliSumOp
import numpy as np

ops = []
h = qubitOp.primitive
h.coeffs = np.real(qubitOp.coeffs)
ops.append(PauliSumOp(h))
print(ops)

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()
options = {"backend" : "ibmq_qasm_simulator"}

from qiskit.circuit.library import RealAmplitudes
circuit = RealAmplitudes(num_qubits=num_qubits, reps=2)

circuit.decompose().draw(fold=100)

from qiskit_ibm_runtime import Estimator
from qiskit.algorithms.optimizers import SPSA

history = {"nfevs": [], "points": [], "fvals": [], "updates": [], "accepted": []}

with Estimator(
    circuits=circuit,
    observables=ops,
    options=options,
    service=service
) as estimator:
    
    def evaluate_expectation(x):
        x = list(x)
        results = estimator(circuits=[0], observables=[0], parameter_values=[x]).values[0]
        return np.real(results)

    def callback(nfev, point, fval, update, accepted):
        print('expectation value: {}'.format(fval))
        history["nfevs"].append(nfev)
        history["points"].append(point)
        history["fvals"].append(fval)
        history["updates"].append(update)
        history["accepted"].append(accepted)    

    np.random.seed(1)
    initial_point = np.random.uniform(-np.pi, np.pi, len(circuit.parameters))
    optimizer = SPSA(30, callback=callback)
    result = optimizer.minimize(fun=evaluate_expectation, x0=initial_point)


print('optimized circuit parameters:\n{}'.format(result.x))

from qiskit_ibm_runtime import Sampler

circut = circuit.measure_all()
with Sampler(
    circuits=circuit,
    options=options,
    service=service
) as sampler:
    sampler_result = sampler(circuits=[0], parameter_values=[result.x], shots=8000)
    print(sampler_result)

bitstrings = sampler_result.quasi_dists[0] # get the dictionary of bitstring counts

import matplotlib.pyplot as plt

fig, axs = plt.subplots(1, 2, figsize=(14, 5))
axs[0].plot(history['fvals'])
axs[1].bar(list(bitstrings.keys()), list(bitstrings.values()))
axs[0].set_xlabel("Energy Evaluation Number")
axs[0].set_ylabel("Energy")
axs[1].set_xlabel("Objective Function Value")
axs[1].set_ylabel("Probability")

samples = []
for x, probability in bitstrings.items():

    x = list(map(int, x)) # convert bitstring to list
    status = qp.is_feasible(x)
    fval = qp.objective.evaluate(x)

    samples.append((x, fval, probability, status))

sorted_samples = sorted(
    samples,
    key=lambda v: (1-v[3], qp.objective.sense.value * v[1]) # feasible: 0; infeasible: 1
)

print('\nsorted')
for sample in sorted_samples:
    print(sample)

best = sorted_samples[0]
print('optimal value: {}'.format(best[0]))
print('optimal function value: {}'.format(best[1]))
print('probability of optimal value: {}'.format(best[2]))

import qiskit_ibm_runtime
qiskit_ibm_runtime.version.get_version_info()

from qiskit.tools.jupyter import *

%qiskit_version_table
%qiskit_copyright