# Path: openfermioncirq/experiments/hfvqe/objective.py
# def generate_hamiltonian(one_body_integrals: np.ndarray,
#                          two_body_integrals: np.ndarray,
#                          constant: float,
#                          EQ_TOLERANCE: Optional[float]=1.0E-12):
#     n_qubits = 2 * one_body_integrals.shape[0]
#     # Initialize Hamiltonian coefficients.
#     one_body_coefficients = np.zeros((n_qubits, n_qubits))
#     two_body_coefficients = np.zeros((n_qubits, n_qubits,
#                                       n_qubits, n_qubits))
#     # Loop through integrals.
#     for p in range(n_qubits // 2):
#         for q in range(n_qubits // 2):
# 
#             # Populate 1-body coefficients. Require p and q have same spin.
#             one_body_coefficients[2 * p, 2 * q] = one_body_integrals[
#                 p, q]
#             one_body_coefficients[2 * p + 1, 2 *
#                                   q + 1] = one_body_integrals[p, q]
#             # Continue looping to prepare 2-body coefficients.
#             for r in range(n_qubits // 2):
#                 for s in range(n_qubits // 2):
#                     # Mixed spin
#                     two_body_coefficients[2 * p, 2 * q + 1,
#                                           2 * r + 1, 2 * s] = (
#                             two_body_integrals[p, q, r, s] / 2.)
#                     two_body_coefficients[2 * p + 1, 2 * q,
#                                           2 * r, 2 * s + 1] = (
#                             two_body_integrals[p, q, r, s] / 2.)
# 
#                     # Same spin
#                     two_body_coefficients[2 * p, 2 * q,
#                                           2 * r, 2 * s] = (
#                             two_body_integrals[p, q, r, s] / 2.)
#                     two_body_coefficients[2 * p + 1, 2 * q + 1,
#                                           2 * r + 1, 2 * s + 1] = (
#                             two_body_integrals[p, q, r, s] / 2.)
# 
#     # Truncate.
#     one_body_coefficients[
#         np.absolute(one_body_coefficients) < EQ_TOLERANCE] = 0.
#     two_body_coefficients[
#         np.absolute(two_body_coefficients) < EQ_TOLERANCE] = 0.
# 
#     # Cast to InteractionOperator class and return.
#     molecular_hamiltonian = InteractionOperator(
#         constant, one_body_coefficients, two_body_coefficients)
#     return molecular_hamiltonian


from typing import Dict, List, Optional
from openfermion.utils import wedge
from openfermion.ops import InteractionRDM
from openfermioncirq.experiments.hfvqe.objective import generate_hamiltonian
import numpy as np
import openfermioncirq.experiments.hfvqe.util as ccu
        opdm_terms = np.random.multivariate_normal(mean=odd_means,
                                                   cov=var_dict['xy_odd'][circuit_idx])
        for idx, (pp0, pp1) in enumerate(o_real_pairs):
            new_opdm[pp0, pp1] = opdm_terms[idx]
            new_opdm[pp1, pp0] = opdm_terms[idx]

        if circuit_idx == 0:
            # resample_diagonal_terms
            opdm_diagonal = np.random.multivariate_normal(
                mean=np.diagonal(opdm), cov=var_dict['z'][circuit_idx])

            # because fill_diagonal documentat seems out of date.
            new_opdm[np.diag_indices_from(new_opdm)] = opdm_diagonal

    return new_opdm


def energy_from_opdm(opdm, constant, one_body_tensor, two_body_tensor):
    """
    Evaluate the energy of an opdm assuming the 2-RDM is opdm ^ opdm

    :param opdm: single spin-component of the full spin-orbital opdm.
    :param constant: constant shift to the Hamiltonian. Commonly this is the
                     nuclear repulsion energy.
    :param one_body_tensor: spatial one-body integrals
    :param two_body_tensor: spatial two-body integrals
    :return:
    """
    spin_opdm = np.kron(opdm, np.eye(2))
    spin_tpdm = 2 * wedge(spin_opdm, spin_opdm, (1, 1), (1, 1))

# Path: openfermioncirq/experiments/hfvqe/objective.py
# def generate_hamiltonian(one_body_integrals: np.ndarray,
#                          two_body_integrals: np.ndarray,
#                          constant: float,
#                          EQ_TOLERANCE: Optional[float]=1.0E-12):
#     n_qubits = 2 * one_body_integrals.shape[0]
#     # Initialize Hamiltonian coefficients.
#     one_body_coefficients = np.zeros((n_qubits, n_qubits))
#     two_body_coefficients = np.zeros((n_qubits, n_qubits,
#                                       n_qubits, n_qubits))
#     # Loop through integrals.
#     for p in range(n_qubits // 2):
#         for q in range(n_qubits // 2):
# 
#             # Populate 1-body coefficients. Require p and q have same spin.
#             one_body_coefficients[2 * p, 2 * q] = one_body_integrals[
#                 p, q]
#             one_body_coefficients[2 * p + 1, 2 *
#                                   q + 1] = one_body_integrals[p, q]
#             # Continue looping to prepare 2-body coefficients.
#             for r in range(n_qubits // 2):
#                 for s in range(n_qubits // 2):
#                     # Mixed spin
#                     two_body_coefficients[2 * p, 2 * q + 1,
#                                           2 * r + 1, 2 * s] = (
#                             two_body_integrals[p, q, r, s] / 2.)
#                     two_body_coefficients[2 * p + 1, 2 * q,
#                                           2 * r, 2 * s + 1] = (
#                             two_body_integrals[p, q, r, s] / 2.)
# 
#                     # Same spin
#                     two_body_coefficients[2 * p, 2 * q,
#                                           2 * r, 2 * s] = (
#                             two_body_integrals[p, q, r, s] / 2.)
#                     two_body_coefficients[2 * p + 1, 2 * q + 1,
#                                           2 * r + 1, 2 * s + 1] = (
#                             two_body_integrals[p, q, r, s] / 2.)
# 
#     # Truncate.
#     one_body_coefficients[
#         np.absolute(one_body_coefficients) < EQ_TOLERANCE] = 0.
#     two_body_coefficients[
#         np.absolute(two_body_coefficients) < EQ_TOLERANCE] = 0.
# 
#     # Cast to InteractionOperator class and return.
#     molecular_hamiltonian = InteractionOperator(
#         constant, one_body_coefficients, two_body_coefficients)
#     return molecular_hamiltonian


# Path: openfermioncirq/experiments/hfvqe/analysis.py
from typing import Dict, List, Optional
from openfermion.utils import wedge
from openfermion.ops import InteractionRDM
from openfermioncirq.experiments.hfvqe.objective import generate_hamiltonian
import numpy as np
import openfermioncirq.experiments.hfvqe.util as ccu

        opdm_terms = np.random.multivariate_normal(mean=odd_means,
                                                   cov=var_dict['xy_odd'][circuit_idx])
        for idx, (pp0, pp1) in enumerate(o_real_pairs):
            new_opdm[pp0, pp1] = opdm_terms[idx]
            new_opdm[pp1, pp0] = opdm_terms[idx]

        if circuit_idx == 0:
            # resample_diagonal_terms
            opdm_diagonal = np.random.multivariate_normal(
                mean=np.diagonal(opdm), cov=var_dict['z'][circuit_idx])

            # because fill_diagonal documentat seems out of date.
            new_opdm[np.diag_indices_from(new_opdm)] = opdm_diagonal

    return new_opdm


def energy_from_opdm(opdm, constant, one_body_tensor, two_body_tensor):
    """
    Evaluate the energy of an opdm assuming the 2-RDM is opdm ^ opdm

    :param opdm: single spin-component of the full spin-orbital opdm.
    :param constant: constant shift to the Hamiltonian. Commonly this is the
                     nuclear repulsion energy.
    :param one_body_tensor: spatial one-body integrals
    :param two_body_tensor: spatial two-body integrals
    :return:
    """
    spin_opdm = np.kron(opdm, np.eye(2))
    spin_tpdm = 2 * wedge(spin_opdm, spin_opdm, (1, 1), (1, 1))

    molecular_hamiltonian = generate_hamiltonian(constant=constant,