How to use the openfermion.ops.QuadOperator function in openfermion

"""
    file_path = get_file_path(file_name, data_directory)

    if plain_text:
        with open(file_path, 'r') as f:
            data = f.read()
            operator_type, operator_terms = data.split(":\n")

        if operator_type == 'FermionOperator':
            operator = FermionOperator(operator_terms)
        elif operator_type == 'BosonOperator':
            operator = BosonOperator(operator_terms)
        elif operator_type == 'QubitOperator':
            operator = QubitOperator(operator_terms)
        elif operator_type == 'QuadOperator':
            operator = QuadOperator(operator_terms)
        else:
            raise TypeError('Operator of invalid type.')
    else:
        with open(file_path, 'rb') as f:
            data = marshal.load(f)
            operator_type = data[0]
            operator_terms = data[1]

        if operator_type == 'FermionOperator':
            operator = FermionOperator()
            for term in operator_terms:
                operator += FermionOperator(term, operator_terms[term])
        elif operator_type == 'BosonOperator':
            operator = BosonOperator()
            for term in operator_terms:
                operator += BosonOperator(term, operator_terms[term])

import numpy as np

from matplotlib import pyplot as plt

from openfermion.ops import QuadOperator
from openfermion.utils import commutator, normal_ordered

import strawberryfields as sf
from strawberryfields.ops import *
from sfopenboson.ops import GaussianPropagation

# define the Hamiltonian
H = QuadOperator('q0 q0', 0.5) + QuadOperator('p0 p0', 0.5) - QuadOperator('q0', 2)

# create the engine
eng = sf.Engine("gaussian")

# set the time-steps
t_vals = np.arange(0, 6, 0.1)
results = np.zeros([2, len(t_vals)])

# evalutate the circuit at each time-step
for step, t in enumerate(t_vals):
    prog = sf.Program(1)

    with prog.context as q:
        Xgate(1) | q[0]
        Zgate(0.5) | q[0]
        GaussianPropagation(H, t) | q

The time evolution unitary associated with the quadratic phase is

    .. math::
        P(s) = \exp\left(i  \frac{s}{2 \hbar} \x^2\right)

    Therefore, :math:`U=e^{-iHt/\hbar}` where
    :math:`H = -\x^2/2` and :math:`t=s`.

    Args:
        s (float): the quadratic phase parameter
        mode (int): the qumode on which the operation acts
    Returns:
        tuple (QuadOperator, t): tuple containing the Hamiltonian
        representing the operation and the propagation time
    """
    return -QuadOperator('q{} q{}'.format(mode, mode))/2, s

.. math::
        CX(s) = \exp\left( -i \frac{s}{\hbar}\x_0\otimes \p_1 \right)

    Therefore, :math:`U=e^{-iHt/\hbar}` where
    :math:`H =\x_0\otimes \p_1` and :math:`t=s`.

    Args:
        s (float): the controlled addition parameter
        mode1 (int): the first qumode :math:`\a_0` on which the operation acts
        mode2 (int): the second qumode :math:`\a_1` on which the operation acts
    Returns:
        tuple (QuadOperator, t): tuple containing the Hamiltonian
        representing the operation and the propagation time
    """
    return QuadOperator('q{} p{}'.format(mode1, mode2)), s

"or use this operator inside an engine context.")
            else:
                self.hbar = hbar

        if not is_hermitian(operator):
            raise ValueError("Hamiltonian must be Hermitian.")

        if (not isinstance(k, int)) or k <= 0:
            raise ValueError("Argument k must be a postive integer.")

        if mode == 'local':
            boson_operator = prune_unused_indices(operator)
        elif mode == 'global':
            boson_operator = operator

        if isinstance(boson_operator, QuadOperator):
            boson_operator = get_boson_operator(boson_operator, hbar=self.hbar)

        self.layer = trotter_layer(boson_operator, t, k)
        self.num_layers = k

        num_modes = max([op[0] for term in operator.terms for op in term])+1

        if mode == 'local':
            self.ns = num_modes
        elif mode == 'global':
            # pylint: disable=protected-access
            self.ns = _Engine._current_context.num_subsystems

a truncation value needs to be provide so that a sparse matrix
    of finite size can be returned.

    Args:
        operator: One of either BosonOperator or QuadOperator.
        trunc (int): The size at which the Fock space should be truncated
            when returning the matrix representing the ladder operator.
        hbar (float): the value of hbar to use in the definition of the
            canonical commutation relation [q_i, p_j] = \delta_{ij} i hbar.
            This only applies if calcualating the sparse representation of
            a quadrature operator.

    Returns:
        The corresponding Scipy sparse matrix of size [trunc, trunc].
    """
    if isinstance(operator, QuadOperator):
        from openfermion.transforms._conversion import get_boson_operator
        boson_operator = get_boson_operator(operator, hbar)
    elif isinstance(operator, BosonOperator):
        boson_operator = operator
    else:
        raise ValueError("Only BosonOperator and QuadOperator are supported.")

    if trunc < 1 or not isinstance(trunc, int):
        raise ValueError("Fock space truncation must be a positive integer.")

    # count the number of modes
    n_modes = 0
    for term in boson_operator.terms:
        for ladder_operator in term:
            if ladder_operator[0] + 1 > n_modes:
                n_modes = ladder_operator[0] + 1

else:
                op = term
                pwr = 1

            mode = int(op[1:])

            if mode not in poly:
                poly[mode] = [0, 0]

            if op[0] == 'q':
                poly[mode][0] += pwr
            elif op[0] == 'p':
                poly[mode][1] += pwr

        # replace {q_i^m p_i^n} -> S({q_i^m p_i^n})
        operator = QuadOperator('')
        for mode, (m, n) in poly.items():
            qtmp = QuadOperator()
            qtmp.terms = mccoy(mode, 'q', 'p', m, n)
            operator *= qtmp
    else:
        operator = QuadOperator.zero()

    return operator

Raises:
        OperatorUtilsError: Not saved, file already exists.
        TypeError: Operator of invalid type.
    """
    file_path = get_file_path(file_name, data_directory)

    if os.path.isfile(file_path) and not allow_overwrite:
        raise OperatorUtilsError("Not saved, file already exists.")

    if isinstance(operator, FermionOperator):
        operator_type = "FermionOperator"
    elif isinstance(operator, BosonOperator):
        operator_type = "BosonOperator"
    elif isinstance(operator, QubitOperator):
        operator_type = "QubitOperator"
    elif isinstance(operator, QuadOperator):
        operator_type = "QuadOperator"
    elif isinstance(operator, (InteractionOperator, InteractionRDM)):
        raise NotImplementedError('Not yet implemented for '
                                  'InteractionOperator or InteractionRDM.')
    else:
        raise TypeError('Operator of invalid type.')

    if plain_text:
        with open(file_path, 'w') as f:
            f.write(operator_type + ":\n" + str(operator))
    else:
        tm = operator.terms
        with open(file_path, 'wb') as f:
            marshal.dump((operator_type,
                          dict(zip(tm.keys(), map(complex, tm.values())))), f)

def get_quad_operator(operator, hbar=1.):
    """Convert to QuadOperator.

    Args:
        operator: BosonOperator.
        hbar (float): the value of hbar used in the definition
            of the commutator [q_i, p_j] = i hbar delta_ij.
            By default hbar=1.

    Returns:
        quad_operator: An instance of the QuadOperator class.
    """
    quad_operator = QuadOperator()

    if isinstance(operator, BosonOperator):
        for term, coefficient in operator.terms.items():
            tmp = QuadOperator('', coefficient)
            for i, d in term:
                tmp *= (1./numpy.sqrt(2.*hbar)) \
                    * (QuadOperator(((i, 'q')))
                        + QuadOperator(((i, 'p')), 1j*(-1)**d))
            quad_operator += tmp

    else:
        raise TypeError("Only BosonOperator is currently "
                        "supported for get_quad_operator.")

    return quad_operator