How to use the openfermion.utils.count_qubits function in openfermion

quadratic_hamiltonian: An instance of the QuadraticHamiltonian class.

    Raises:
        TypeError: Input must be a FermionOperator.
        TypeError: FermionOperator does not map to QuadraticHamiltonian.

    Warning:
        Even assuming that each creation or annihilation operator appears
        at most a constant number of times in the original operator, the
        runtime of this method is exponential in the number of qubits.
    """
    if not isinstance(fermion_operator, FermionOperator):
        raise TypeError('Input must be a FermionOperator.')

    if n_qubits is None:
        n_qubits = count_qubits(fermion_operator)
    if n_qubits < count_qubits(fermion_operator):
        raise ValueError('Invalid number of qubits specified.')

    # Normal order the terms and initialize.
    fermion_operator = normal_ordered(fermion_operator)
    constant = 0.
    combined_hermitian_part = numpy.zeros((n_qubits, n_qubits), complex)
    antisymmetric_part = numpy.zeros((n_qubits, n_qubits), complex)

    # Loop through terms and assign to matrix.
    for term in fermion_operator.terms:
        coefficient = fermion_operator.terms[term]
        # Ignore this term if the coefficient is zero
        if abs(coefficient) < EQ_TOLERANCE:
            continue

original_ordering (list): The initial Jordan-Wigner canonical order.

    Returns:
        A 3-tuple of terms from the Hubbard Hamiltonian in order of
        simulation, the indices they act on, and whether they are hopping
        operators (both also in the same order).

    Notes:
        Assumes that the Hubbard model has spin and is on a 2D square
        aperiodic lattice. Uses the "stagger"-based Trotter step for the
        Hubbard model detailed in Kivlichan et al., "Quantum Simulation
        of Electronic Structure with Linear Depth and Connectivity",
        arxiv:1711.04789.
    """
    hamiltonian = normal_ordered(hubbard_hamiltonian)
    n_qubits = count_qubits(hamiltonian)
    side_length = int(numpy.sqrt(n_qubits / 2.0))

    ordered_terms = []
    ordered_indices = []
    ordered_is_hopping_operator = []

    original_ordering = list(range(n_qubits))
    for i in range(0, n_qubits - side_length, 2 * side_length):
        for j in range(2 * bool(i % (4 * side_length)), 2 * side_length, 4):
            original_ordering[i+j], original_ordering[i+j+1] = (
                original_ordering[i+j+1], original_ordering[i+j])

    input_ordering = list(original_ordering)

    # Follow odd-even transposition sort. In alternating steps, swap each even
    # qubits with the odd qubit to its right, and in the next step swap each

def bravyi_kitaev_fast_edge_matrix(iop, n_qubits=None):
    """
    Use InteractionOperator to construct edge matrix required for the algorithm

    Edge matrix contains the information about the edges between vertices.
    Edge matrix is required to build the operators in bravyi_kitaev_fast model.

    Args:
        iop (Interaction Operator):

    Returns:
        edge_matrix (Numpy array): A square numpy array containing information
            about the edges present in the model.
    """
    n_qubits = count_qubits(iop)
    edge_matrix = 1j*numpy.zeros((n_qubits, n_qubits))
    # Loop through all indices.
    for p in range(n_qubits):
        for q in range(n_qubits):

            # Handle one-body terms.
            coefficient = complex(iop[(p, 1), (q, 0)])
            if coefficient and p >= q:
                edge_matrix[p, q] = bool(complex(iop[(p, 1), (q, 0)]))

            # Keep looping for the two-body terms.
            for r in range(n_qubits):
                for s in range(n_qubits):
                    coefficient2 = complex(iop[(p, 1), (q, 1), (r, 0), (s, 0)])

                    # Skip zero terms.

n_qubits: the number of qubits term acts on. If not set, defaults
                to the maximum qubit number acted on by term.

    Returns:
        transformed_term: An instance of the FermionOperator class.

    Raises:
        TypeError: Input must be a QubitOperator.
        TypeError: Invalid number of qubits specified.
        TypeError: Pauli operators must be X, Y or Z.
    """
    if not isinstance(qubit_operator, QubitOperator):
        raise TypeError('Input must be a QubitOperator.')
    if n_qubits is None:
        n_qubits = count_qubits(qubit_operator)
    if n_qubits < count_qubits(qubit_operator):
        raise ValueError('Invalid number of qubits specified.')

    # Loop through terms.
    transformed_operator = FermionOperator()
    for term in qubit_operator.terms:
        transformed_term = FermionOperator(())
        if term:
            working_term = QubitOperator(term)
            pauli_operator = term[-1]
            while pauli_operator is not None:

                # Handle Pauli Z.
                if pauli_operator[1] == 'Z':
                    transformed_pauli = FermionOperator(
                        ()) + number_operator(n_qubits, pauli_operator[0], -2.)

def jordan_wigner_diagonal_coulomb_hamiltonian(operator):
    n_qubits = count_qubits(operator)
    qubit_operator = QubitOperator((), operator.constant)

    # Transform diagonal one-body terms
    for p in range(n_qubits):
        coefficient = operator.one_body[p, p] + operator.two_body[p, p]
        qubit_operator += QubitOperator(((p, 'Z'),), -.5 * coefficient)
        qubit_operator += QubitOperator((), .5 * coefficient)

    # Transform other one-body terms and two-body terms
    for p, q in itertools.combinations(range(n_qubits), 2):
        # One-body
        real_part = numpy.real(operator.one_body[p, q])
        imag_part = numpy.imag(operator.one_body[p, q])
        parity_string = [(i, 'Z') for i in range(p + 1, q)]
        qubit_operator += QubitOperator(
                [(p, 'X')] + parity_string + [(q, 'X')], .5 * real_part)

def _bravyi_kitaev_fermion_operator(operator, n_qubits):
    # Compute the number of qubits.
    N = count_qubits(operator)
    if n_qubits is None:
        n_qubits = N
    if n_qubits < N:
        raise ValueError('Invalid number of qubits specified.')

    # Compute transformed operator.
    transformed_terms = (
        _transform_operator_term(term=term,
                                 coefficient=operator.terms[term],
                                 n_qubits=n_qubits)
        for term in operator.terms
    )
    return inline_sum(summands=transformed_terms, seed=QubitOperator())

A FermionOperator to transform.
        n_qubits (int|None):
            Can force the number of qubits in the resulting operator above the
            number that appear in the input operator.

    Returns:
        transformed_operator: An instance of the QubitOperator class.

    Raises:
        ValueError: Invalid number of qubits specified.
    """
    # Compute the number of qubits.
    from openfermion.utils import count_qubits
    if n_qubits is None:
        n_qubits = count_qubits(operator)
    if n_qubits < count_qubits(operator):
        raise ValueError('Invalid number of qubits specified.')

    # Build the Fenwick tree.
    fenwick_tree = FenwickTree(n_qubits)

    # Compute transformed operator.
    transformed_terms = (
        _transform_operator_term(term=term,
                                 coefficient=operator.terms[term],
                                 fenwick_tree=fenwick_tree)
        for term in operator.terms
    )
    return inline_sum(summands=transformed_terms, seed=QubitOperator())

interaction_operator: An instance of the InteractionOperator class.

    Raises:
        TypeError: Input must be a FermionOperator.
        TypeError: FermionOperator does not map to InteractionOperator.

    Warning:
        Even assuming that each creation or annihilation operator appears
        at most a constant number of times in the original operator, the
        runtime of this method is exponential in the number of qubits.
    """
    if not isinstance(fermion_operator, FermionOperator):
        raise TypeError('Input must be a FermionOperator.')

    if n_qubits is None:
        n_qubits = count_qubits(fermion_operator)
    if n_qubits < count_qubits(fermion_operator):
        raise ValueError('Invalid number of qubits specified.')

    # Normal order the terms and initialize.
    fermion_operator = normal_ordered(fermion_operator)
    constant = 0.
    one_body = numpy.zeros((n_qubits, n_qubits), complex)
    two_body = numpy.zeros((n_qubits, n_qubits,
                            n_qubits, n_qubits), complex)

    # Loop through terms and assign to matrix.
    for term in fermion_operator.terms:
        coefficient = fermion_operator.terms[term]
        # Ignore this term if the coefficient is zero
        if abs(coefficient) < EQ_TOLERANCE:
            continue

def apply_constraints(operator, n_fermions):
    """Function to use linear programming to apply constraints.

    Args:
        operator(FermionOperator): FermionOperator with only 1- and 2-body
            terms that we wish to vectorize.
        n_fermions(int): The number of particles in the simulation.

    Returns:
        modified_operator(FermionOperator): The operator with reduced norm
            that has been modified with equality constraints.
    """
    # Get constraint matrix.
    n_orbitals = count_qubits(operator)
    constraints = constraint_matrix(n_orbitals, n_fermions)
    n_constraints, n_terms = constraints.get_shape()

    # Get vectorized operator.
    vectorized_operator = operator_to_vector(operator)
    initial_bound = numpy.sum(numpy.absolute(vectorized_operator[1::])) ** 2
    print('Initial bound on measurements is %f.' % initial_bound)

    # Get linear programming coefficient vector.
    n_variables = n_constraints + n_terms
    lp_vector = numpy.zeros(n_variables, float)
    lp_vector[-n_terms:] = 1.

    # Get linear programming constraint matrix.
    lp_constraint_matrix = scipy.sparse.dok_matrix((2 * n_terms, n_variables))
    for (i, j), value in constraints.items():

def lowest_single_particle_energy_states(hamiltonian, n_states):
    """Find the lowest single-particle states of the given Hamiltonian.

    Args:
        hamiltonian (FermionOperator)
        n_states (int): Number of lowest energy states to give."""
    # Enumerate the single-particle states.
    n_single_particle_states = count_qubits(hamiltonian)

    # Compute the energies for each of the single-particle states.
    single_particle_energies = numpy.zeros(n_single_particle_states,
                                           dtype=float)
    for i in range(n_single_particle_states):
        single_particle_energies[i] = hamiltonian.terms.get(
            ((i, 1), (i, 0)), 0.0)

    # Find the n_states lowest states.
    occupied_states = single_particle_energies.argsort()[:n_states]

    return list(occupied_states)