import numpy as np
from qiskit_nature.operators.second_quantization import QuadraticHamiltonian

# create Hamiltonian
hermitian_part = np.array(
    [
        [1.0, 2.0, 0.0, 0.0],
        [2.0, 1.0, 2.0, 0.0],
        [0.0, 2.0, 1.0, 2.0],
        [0.0, 0.0, 2.0, 1.0],
    ]
)
antisymmetric_part = np.array(
    [
        [0.0, 3.0, 0.0, 0.0],
        [-3.0, 0.0, 3.0, 0.0],
        [0.0, -3.0, 0.0, 3.0],
        [0.0, 0.0, -3.0, 0.0],
    ]
)
constant = 4.0

hamiltonian = QuadraticHamiltonian(
    hermitian_part=hermitian_part,
    antisymmetric_part=antisymmetric_part,
    constant=constant,
)

# convert it to a FermionicOp and print it
hamiltonian_ferm = hamiltonian.to_fermionic_op()
hamiltonian_ferm.display_format = "sparse"
hamiltonian_ferm.set_truncation(0)
print(hamiltonian_ferm)

# get the transformation matrix W and orbital energies {epsilon_j}
(
    transformation_matrix,
    orbital_energies,
    transformed_constant,
) = hamiltonian.diagonalizing_bogoliubov_transform()

print("Shape of matrix W:")
print(transformation_matrix.shape)
print()
print("Orbital energies:")
print(orbital_energies)
print()
print("Transformed constant:")
print(transformed_constant)

from qiskit_nature.circuit.library import FermionicGaussianState

occupied_orbitals = (0, 2)
eig = np.sum(orbital_energies[list(occupied_orbitals)]) + transformed_constant
print("Eigenvalue:")
print(eig)

circuit = FermionicGaussianState(
    transformation_matrix, occupied_orbitals=occupied_orbitals
)
circuit.draw("mpl")

from qiskit.quantum_info import Statevector
from qiskit_nature.mappers.second_quantization import JordanWignerMapper

# simulate the circuit to get the final state
state = np.array(Statevector(circuit))

# convert the Hamiltonian to a matrix
hamiltonian_jw = JordanWignerMapper().map(hamiltonian_ferm).to_matrix()

# check that the state is an eigenvector with the expected eigenvalue
np.testing.assert_allclose(hamiltonian_jw @ state, eig * state, atol=1e-8)

# create Hamiltonian
hermitian_part = np.array(
    [
        [1.0, 2.0, 0.0, 0.0],
        [2.0, 1.0, 2.0, 0.0],
        [0.0, 2.0, 1.0, 2.0],
        [0.0, 0.0, 2.0, 1.0],
    ]
)
constant = 4.0

hamiltonian = QuadraticHamiltonian(
    hermitian_part=hermitian_part,
    constant=constant,
)

print("Hamiltonian conserves particle number:")
print(hamiltonian.conserves_particle_number())

# get the transformation matrix W and orbital energies {epsilon_j}
(
    transformation_matrix,
    orbital_energies,
    transformed_constant,
) = hamiltonian.diagonalizing_bogoliubov_transform()

print("Shape of matrix W:")
print(transformation_matrix.shape)
print()
print("Orbital energies:")
print(orbital_energies)
print()
print("Transformed constant:")
print(transformed_constant)

from qiskit_nature.circuit.library import SlaterDeterminant

occupied_orbitals = (0, 2)
eig = np.sum(orbital_energies[list(occupied_orbitals)]) + transformed_constant
print("Eigenvalue:")
print(eig)

circuit = SlaterDeterminant(transformation_matrix[list(occupied_orbitals)])
circuit.draw("mpl")

from qiskit.quantum_info import random_hermitian

# construct and display some larger circuits
N = 8
n_particles = 4
hermitian_part = np.array(random_hermitian(N))
hamiltonian = QuadraticHamiltonian(hermitian_part)
transformation_matrix, _, _ = hamiltonian.diagonalizing_bogoliubov_transform()
occupied_orbitals = range(n_particles)
circuit = SlaterDeterminant(transformation_matrix[list(occupied_orbitals)])
circuit.draw("mpl")