OpenFermion: Bad test for Bravyi-Kitaev SuperFast Transform

The file transforms/_bksf_test.py contains the following test:

    def test_bravyi_kitaev_fast_number_excitation_operator(self):
        # using hydrogen Hamiltonian and introducing some number operator terms
        constant = 0
        one_body = numpy.zeros((4, 4))
        one_body[(0, 0)] = .4
        one_body[(1, 1)] = .5
        one_body[(2, 2)] = .6
        one_body[(3, 3)] = .7
        two_body = self.molecular_hamiltonian.two_body_tensor
        # initiating number operator terms for all the possible cases
        two_body[(1, 2, 3, 1)] = 0.1
        two_body[(1, 3, 2, 1)] = 0.1
        two_body[(1, 2, 1, 3)] = 0.15
        two_body[(3, 1, 2, 1)] = 0.15
        two_body[(0, 2, 2, 1)] = 0.09
        two_body[(1, 2, 2, 0)] = 0.09
        two_body[(1, 2, 3, 2)] = 0.11
        two_body[(2, 3, 2, 1)] = 0.11
        two_body[(2, 2, 2, 2)] = 0.1
        molecular_hamiltonian = InteractionOperator(constant,
                                                    one_body, two_body)
        # comparing the eigenspectrum of Hamiltonian
        n_qubits = count_qubits(molecular_hamiltonian)
        bravyi_kitaev_fast_H = _bksf.bravyi_kitaev_fast(molecular_hamiltonian)
        jw_H = jordan_wigner(molecular_hamiltonian)
        bravyi_kitaev_fast_H_eig = eigenspectrum(bravyi_kitaev_fast_H)
        jw_H_eig = eigenspectrum(jw_H)
        bravyi_kitaev_fast_H_eig = bravyi_kitaev_fast_H_eig.round(5)
        jw_H_eig = jw_H_eig.round(5)
        evensector_H = 0
        for i in range(numpy.size(jw_H_eig)):
            if bool(numpy.size(numpy.where(jw_H_eig[i] ==
                                           bravyi_kitaev_fast_H_eig))):
                evensector_H += 1

        # comparing eigenspectrum of number operator
        bravyi_kitaev_fast_n = _bksf.number_operator(molecular_hamiltonian)
        jw_n = QubitOperator()
        n_qubits = count_qubits(molecular_hamiltonian)
        for i in range(n_qubits):
            jw_n += jordan_wigner_one_body(i, i)
        jw_eig_spec = eigenspectrum(jw_n)
        bravyi_kitaev_fast_eig_spec = eigenspectrum(bravyi_kitaev_fast_n)
        evensector_n = 0
        for i in range(numpy.size(jw_eig_spec)):
            if bool(numpy.size(numpy.where(jw_eig_spec[i] ==
                                           bravyi_kitaev_fast_eig_spec))):
                evensector_n += 1
        self.assertEqual(evensector_H, 2**(n_qubits - 1))
        self.assertEqual(evensector_n, 2**(n_qubits - 1))

The line

jw_H = jordan_wigner(molecular_hamiltonian)

depends on the Jordan-Wigner routine for InteractionOperators, which is incorrect, as pointed out in #281 . Indeed, if we change this line to use the standard Jordan-Wigner,

jw_H = jordan_wigner(get_fermion_operator(molecular_hamiltonian))

then the test fails. At the end, evensector_H has value 2 instead of 8. This is the same reason that my PR #286 cannot pass Travis.

The routines get_fermion_operator and jordan_wigner applied to FermionOperators are both well-tested and have simple source code so I have high confidence in them.

I was hoping somebody more familiar with the BKSF code could help us figure out what is going on here. Perhaps @kanavsetia?