# Path: openfermioncirq/experiments/hfvqe/circuits.py
# def rhf_params_to_matrix(parameters: np.ndarray, num_qubits: int,
#                          occ: Optional[Union[None, List[int]]]=None,
#                          virt: Optional[Union[None, List[int]]]=None):
#     """
#     For restricted Hartree-Fock we have nocc * nvirt parameters.  These are provided
#     as a list that is ordered by (virtuals) \times (occupied) where
#     occupied is a set of indices corresponding to the occupied oribitals w.r.t the
#     Lowdin basis and virtuals is a set of indices of the virutal orbitals w.r.t the
#     Lowdin basis.  For example, for H4 we have 2 orbitals occupied and 2 virtuals
# 
#     occupied = [0, 1]  virtuals = [2, 3]
# 
#     parameters = [(v_{0}, o_{0}), (v_{0}, o_{1}), (v_{1}, o_{0}), (v_{1}, o_{1})]
#                = [(2, 0), (2, 1), (3, 0), (3, 1)]
# 
#     You can think of the tuples of elements of the upper right triangle of the
#     antihermitian matrix that specifies the c_{b, i} coefficients.
# 
#     coefficient matrix
#     [[ c_{0, 0}, -c_{1, 0}, -c_{2, 0}, -c_{3, 0}],
#      [ c_{1, 0},  c_{1, 1}, -c_{2, 1}, -c_{3, 1}],
#      [ c_{2, 0},  c_{2, 1},  c_{2, 2}, -c_{3, 2}],
#      [ c_{3, 0},  c_{3, 1},  c_{3, 2},  c_{3, 3}]]
# 
#     Since we are working with only non-redundant operators we know c_{i, i} = 0
#     and any c_{i, j} where i and j are both in occupied or both in virtual = 0.
#     """
#     if occ is None:
#         occ = range(num_qubits//2)
#     if virt is None:
#         virt = range(num_qubits // 2, num_qubits)
# 
#     # check that parameters are a real array
#     if not np.allclose(parameters.imag, 0):
#         raise ValueError("parameters input must be real valued")
# 
#     kappa = np.zeros((len(occ) + len(virt), len(occ) + len(virt)))
#     for idx, (v, o) in enumerate(product(virt, occ)):
#         kappa[v, o] = parameters[idx].real
#         kappa[o, v] = -parameters[idx].real
#     return kappa


from typing import Optional
from itertools import product
from openfermion.ops import InteractionOperator, InteractionRDM
from openfermion.utils import wedge
from openfermion.transforms import get_fermion_operator
from openfermioncirq.experiments.hfvqe.circuits import rhf_params_to_matrix
import numpy as np

        :param opdm:
        :return:
        """
        rdms = self.rdms_from_opdm_aa(opdm_aa)
        return rdms.expectation(self.hamiltonian).real

    def global_gradient_opdm(self, params: np.ndarray,
                              alpha_opdm: np.ndarray):
        opdm = np.zeros((self.num_qubits, self.num_qubits),
                           dtype=np.complex128)
        opdm[::2, ::2] = alpha_opdm
        opdm[1::2, 1::2] = alpha_opdm
        tpdm = 2 * wedge(opdm, opdm, (1, 1), (1, 1))

        # now go through and generate all the necessary Z, Y, Y_kl matrices
        kappa_matrix = rhf_params_to_matrix(params,
                                            len(self.occ) + len(self.virt),
                                            self.occ, self.virt)
        kappa_matrix_full = np.kron(kappa_matrix, np.eye(2))
        w_full, v_full = np.linalg.eigh(-1j * kappa_matrix_full)  # so that kappa = i U lambda U^
        eigs_scaled_full = get_matrix_of_eigs(w_full)

        grad = np.zeros(self.nocc * self.nvirt, dtype=np.complex128)
        kdelta = np.eye(self.num_qubits)

        # NOW GENERATE ALL TERMS ASSOCIATED WITH THE GRADIENT!!!!!!
        for p in range(self.nocc * self.nvirt):
            grad_params = np.zeros_like(params)
            grad_params[p] = 1

# Path: openfermioncirq/experiments/hfvqe/circuits.py
# def rhf_params_to_matrix(parameters: np.ndarray, num_qubits: int,
#                          occ: Optional[Union[None, List[int]]]=None,
#                          virt: Optional[Union[None, List[int]]]=None):
#     """
#     For restricted Hartree-Fock we have nocc * nvirt parameters.  These are provided
#     as a list that is ordered by (virtuals) \times (occupied) where
#     occupied is a set of indices corresponding to the occupied oribitals w.r.t the
#     Lowdin basis and virtuals is a set of indices of the virutal orbitals w.r.t the
#     Lowdin basis.  For example, for H4 we have 2 orbitals occupied and 2 virtuals
# 
#     occupied = [0, 1]  virtuals = [2, 3]
# 
#     parameters = [(v_{0}, o_{0}), (v_{0}, o_{1}), (v_{1}, o_{0}), (v_{1}, o_{1})]
#                = [(2, 0), (2, 1), (3, 0), (3, 1)]
# 
#     You can think of the tuples of elements of the upper right triangle of the
#     antihermitian matrix that specifies the c_{b, i} coefficients.
# 
#     coefficient matrix
#     [[ c_{0, 0}, -c_{1, 0}, -c_{2, 0}, -c_{3, 0}],
#      [ c_{1, 0},  c_{1, 1}, -c_{2, 1}, -c_{3, 1}],
#      [ c_{2, 0},  c_{2, 1},  c_{2, 2}, -c_{3, 2}],
#      [ c_{3, 0},  c_{3, 1},  c_{3, 2},  c_{3, 3}]]
# 
#     Since we are working with only non-redundant operators we know c_{i, i} = 0
#     and any c_{i, j} where i and j are both in occupied or both in virtual = 0.
#     """
#     if occ is None:
#         occ = range(num_qubits//2)
#     if virt is None:
#         virt = range(num_qubits // 2, num_qubits)
# 
#     # check that parameters are a real array
#     if not np.allclose(parameters.imag, 0):
#         raise ValueError("parameters input must be real valued")
# 
#     kappa = np.zeros((len(occ) + len(virt), len(occ) + len(virt)))
#     for idx, (v, o) in enumerate(product(virt, occ)):
#         kappa[v, o] = parameters[idx].real
#         kappa[o, v] = -parameters[idx].real
#     return kappa


# Path: openfermioncirq/experiments/hfvqe/objective.py
from typing import Optional
from itertools import product
from openfermion.ops import InteractionOperator, InteractionRDM
from openfermion.utils import wedge
from openfermion.transforms import get_fermion_operator
from openfermioncirq.experiments.hfvqe.circuits import rhf_params_to_matrix
import numpy as np


        :param opdm:
        :return:
        """
        rdms = self.rdms_from_opdm_aa(opdm_aa)
        return rdms.expectation(self.hamiltonian).real

    def global_gradient_opdm(self, params: np.ndarray,
                              alpha_opdm: np.ndarray):
        opdm = np.zeros((self.num_qubits, self.num_qubits),
                           dtype=np.complex128)
        opdm[::2, ::2] = alpha_opdm
        opdm[1::2, 1::2] = alpha_opdm
        tpdm = 2 * wedge(opdm, opdm, (1, 1), (1, 1))

        # now go through and generate all the necessary Z, Y, Y_kl matrices
        kappa_matrix = rhf_params_to_matrix(params,
                                            len(self.occ) + len(self.virt),
                                            self.occ, self.virt)
        kappa_matrix_full = np.kron(kappa_matrix, np.eye(2))
        w_full, v_full = np.linalg.eigh(-1j * kappa_matrix_full)  # so that kappa = i U lambda U^
        eigs_scaled_full = get_matrix_of_eigs(w_full)

        grad = np.zeros(self.nocc * self.nvirt, dtype=np.complex128)
        kdelta = np.eye(self.num_qubits)

        # NOW GENERATE ALL TERMS ASSOCIATED WITH THE GRADIENT!!!!!!
        for p in range(self.nocc * self.nvirt):
            grad_params = np.zeros_like(params)
            grad_params[p] = 1

            Y = rhf_params_to_matrix(grad_params,