"""Variational Quantum Eigensolver (VQE) tutorial."""

from IPython.display import clear_output
import matplotlib.pyplot as plt
import numpy as np

from qiskit_nature.units import DistanceUnit
from qiskit_nature.second_q.drivers import PySCFDriver
from qiskit_nature.second_q.mappers import JordanWignerMapper
from qiskit_nature.second_q.circuit.library import HartreeFock
from qiskit_nature.second_q.algorithms import GroundStateEigensolver
from qiskit_nature import settings

from qiskit_algorithms.utils import algorithm_globals
from qiskit.circuit.library import EfficientSU2
from qiskit_algorithms.optimizers import L_BFGS_B
from qiskit_algorithms.minimum_eigensolvers import NumPyMinimumEigensolver, VQE

from qiskit_qulacs.qulacs_estimator import QulacsEstimator
from qiskit_qulacs.qulacs_estimator_gradient import QulacsEstimatorGradient

settings.use_pauli_sum_op = False

algorithm_globals.random_seed = 42

driver = PySCFDriver(
    atom="H 0 0 0; H 0 0 0.735",
    basis="sto3g",
    charge=0,
    spin=0,
    unit=DistanceUnit.ANGSTROM,
)

problem = driver.run()
fermionic_op = problem.hamiltonian.second_q_op()

mapper = JordanWignerMapper()
qubit_op = mapper.map(fermionic_op)
print(qubit_op)

print("Number of qubits:", qubit_op.num_qubits)
print("Number of hamiltonian terms:", len(qubit_op))

# Get reference solution
numpy_solver = NumPyMinimumEigensolver()
calc = GroundStateEigensolver(mapper, numpy_solver)
res_actual = calc.solve(problem)
print(res_actual)

exact_energy = res_actual.eigenvalues[0]
print(f"Expected ground state energy: {exact_energy:.12f} Ha")

# Hartree focks state
num_particles = problem.num_particles
num_spatial_orbitals = problem.num_spatial_orbitals

init_state = HartreeFock(
    num_spatial_orbitals=num_spatial_orbitals,
    num_particles=num_particles,
    qubit_mapper=mapper,
)

# ansatz
ansatz = EfficientSU2(qubit_op.num_qubits, su2_gates=["ry"]).decompose()

# Add the initial state
init_state.barrier()
ansatz = init_state.compose(ansatz)
ansatz.draw("mpl")

intermediate_info = []


def callback(eval_count, parameters, value, std):  # pylint: disable=unused-argument
    """
    A callback that can be registered with the optimizer to store the intermediate
    value and parameters.
    """
    intermediate_info.append(value)
    clear_output(wait=True)
    plt.plot(
        intermediate_info,
        color="purple",
        lw=2,
        label=f"Simulated VQE {np.round(value,4)}",
    )
    plt.ylabel("Energy")
    plt.xlabel("Iterations")

    # Exact ground state energy value
    plt.axhline(
        y=exact_energy,
        color="tab:red",
        ls="--",
        lw=2,
        label="Target: " + str(np.round(exact_energy, 4)),
    )
    plt.legend()
    plt.grid()
    plt.show()

optimizer = L_BFGS_B(maxiter=20)

qulacs_estimator = QulacsEstimator()
qulacs_gradient = QulacsEstimatorGradient()

vqe = VQE(
    qulacs_estimator, ansatz, optimizer, callback=callback, gradient=qulacs_gradient
)
result = vqe.compute_minimum_eigenvalue(operator=qubit_op)

print(result)

def rel_err(target, measured):
    """Compute the relative error."""
    return abs((target - measured) / target)

# Compute the relative error between the expected ground state energy and the VQE's output
rel_error = rel_err(exact_energy, result.eigenvalue)

print(f"Expected ground state energy: {exact_energy:.12f}")
print(f"Computed ground state energy: {result.eigenvalue:.12f}")
print(f"Relative error: {rel_error:.12f}")