!pip install qiskit-algorithms

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService(channel="ibm_quantum")
backend = service.backend("ibmq_qasm_simulator")

from qiskit.quantum_info import SparsePauliOp

# Create an operator by summing together three Pauli operators 
choc_op = SparsePauliOp(['ZII', 'IZI', 'IIZ'])

print(choc_op.to_matrix())


import numpy as np
from qiskit_algorithms.optimizers import SPSA
from qiskit.circuit.library import EfficientSU2
from qiskit_ibm_runtime import Session, Estimator

ansatz = EfficientSU2(3)
ansatz.decompose().draw('mpl')

from qiskit_ibm_runtime import Options

options = Options()
options.resilience_level = 1
options.execution.shots = 1024

with Session(backend=backend, max_time="1h") as session:
    estimator = Estimator(options=options)
    
    # Define a callback to print values in each iteration of the optimizer
    callback = lambda nfev, params, fval, step, acc: print(f'Evaluation {nfev}: {fval}')
    optimizer = SPSA(callback=callback, maxiter=800)

    # Define a cost function, flipping the sign so that the optimizer maximizes the expectation value
    cost_func = lambda params: estimator.run(ansatz, choc_op, parameter_values=params).result().values[0] * -1
    result = optimizer.minimize(cost_func, x0=np.zeros(ansatz.num_parameters))

print(result)

import matplotlib.pyplot as plt

jobs = service.jobs(session_id=session.session_id, limit=None)
expectations = [job.result().values[0] for job in jobs if job.done()][::-1]

plt.plot(range(len(expectations)), expectations, label='backend')
plt.plot(range(len(expectations)), [3.0]*len(expectations), ls='--', label='exact')
plt.xlabel('Iteration', fontsize=16)
plt.ylabel('Expectation', fontsize=16)
plt.legend(loc='lower right')

