Repository: quantumlib/OpenFermion
Pull request #516 (2019-05-15 13:13:32)

Tapering off qubits using stabilizer conditions.
This module tapers off qubits given a set of stabilizers.
Additionally, it reduces the number of terms of a Hamiltonian (operator) using stabilizer conditions.

The module allows the user to:
1) Select which qubit to taper off by a given stabilizer.
2) Keep the length of the Pauli string after reducing terms the same as the starting one.
3) Return which qubits have been tapered off.

The module comes with documentation and test module.

Xavier Bonet-Monroig and Mark Steudtner

googlebot: 
Thanks for your pull request. It looks like this may be your first contribution to a Google open source project (if not, look below for help). Before we can look at your pull request, you'll need to sign a Contributor License Agreement (CLA).

:memo: **Please visit <https://cla.developers.google.com/> to sign.**

Once you've signed (or fixed any issues), please reply here (e.g. `I signed it!`) and we'll verify it.

----

#### What to do if you already signed the CLA

##### Individual signers

*   It's possible we don't have your GitHub username or you're using a different email address on your commit. Check [your existing CLA data](https://cla.developers.google.com/clas) and verify that your [email is set on your git commits](https://help.github.com/articles/setting-your-email-in-git/).

##### Corporate signers

*   Your company has a Point of Contact who decides which employees are authorized to participate. Ask your POC to be added to the group of authorized contributors. If you don't know who your Point of Contact is, direct the Google project maintainer to [go/cla#troubleshoot](http://go/cla#troubleshoot) ([Public version](https://opensource.google.com/docs/cla/#troubleshoot)).
*   The email used to register you as an authorized contributor must be the email used for the Git commit. Check [your existing CLA data](https://cla.developers.google.com/clas) and verify that your [email is set on your git commits](https://help.github.com/articles/setting-your-email-in-git/).
*   The email used to register you as an authorized contributor must also be [attached to your GitHub account](https://github.com/settings/emails).
		

ℹ️ **Googlers: [Go here](https://goto.google.com/prinfo/https%3A%2F%2Fgithub.com%2Fquantumlib%2FOpenFermion%2Fpull%2F516) for more info**.

<!-- need_sender_cla -->
xabomon: I have signed the CLA!
googlebot: So there's good news and bad news.

:thumbsup: The good news is that everyone that needs to sign a CLA (the pull request submitter and all commit authors) have done so.  Everything is all good there.

:confused: The bad news is that it appears that one or more commits were authored or co-authored by someone other than the pull request submitter.  We need to confirm that all authors are ok with their commits being contributed to this project.  Please have them confirm that here in the pull request.

*Note to project maintainer: This is a terminal state, meaning the `cla/google` commit status will not change from this state. It's up to you to confirm consent of all the commit author(s), set the `cla` label to `yes` (if enabled on your project), and then merge this pull request when appropriate.*

ℹ️ **Googlers: [Go here](https://goto.google.com/prinfo/https%3A%2F%2Fgithub.com%2Fquantumlib%2FOpenFermion%2Fpull%2F516) for more info**.

<!-- need_author_consent -->
msteudtner: I'm ok with my comments being contributed.
ncrubin: Hey @xabomon thanks for the PR.  Looks great.  Can you merge with master so we get a fresh CI build and we can check out why this is failing. I suspect something with the init registration.  
ncrubin: @xabomon looks like on of your tests is still attempting to import `reduce_number_of_terms` from transforms and not utils.
xabomon: Hi @ncrubin,

Sorry for this, I just fixed the import in the test module.

Best,
Xavi
xabomon: Hi @ncrubin ,

After yesterday's commit I saw in the CI that one of the test was failing in the python2.7 version because I was defining it as an empty list, after changing the mutable variables to None.

I updated the test module that includes the new possible case, and change the old one. This should now have no problem in the CI.

Let me know if anything else needs to be updated.

Best,
Xavi
ncrubin: I believe the rules for implementing `super()` are different in python2.7.  I can look this weekend to resolve this if you don't solve by then!
xabomon: Hi Nick,

Don’t worry about it! I can’t look at it during the weekend, but I will do it next week.

Thanks,
Xavi
ncrubin: @xabomon bump on this?  it would be great to get this in!
xabomon: Hi Nick,

I'm really sorry for not being active on this lately.

I manage to narrow down that the test fail in python 2.7. I am suspicious that the fail is due to our customize error messages but I have seen code within OpenFermion with customize errors.

Currently I'm building an environment with 2.7 to see if I can reproduce the issue outside the test module.

Please, let me know if you have heard about similar issues before and how to address them.

Best,
Xavi
xabomon: Hi @ncrubin ,
We have fixed the CI issue but now the coverage is failing.

As far as I understand, the coveralls counts how many new lines are used in the test module. We have tested all possible issue in the main functions, this implies that internal functions are indirectly tested. Do we have to also cover internal functions?

Thanks for the assistance.

Xavi
xabomon: I finally managed to fulfill CI and coverage.

Hopefully, we can now merge this PR.

Best,
Xavi
googlebot: :frowning_face: Sorry, but only Googlers may change the label `cla: yes`.
googlebot: :frowning_face: Sorry, but only Googlers may change the label `cla: yes`.
googlebot: :frowning_face: Sorry, but only Googlers may change the label `cla: yes`.

Patch:
diff --git a/src/openfermion/utils/__init__.py b/src/openfermion/utils/__init__.py
index 8f6454c..bf863ef 100644
--- a/src/openfermion/utils/__init__.py
+++ b/src/openfermion/utils/__init__.py
@@ -35,6 +35,9 @@ from ._operator_utils import (chemist_ordered, count_qubits,
                               load_operator, save_operator,
                               group_into_tensor_product_basis_sets)
 
+from ._qubit_tapering_from_stabilizer import (reduce_number_of_terms,
+                                              taper_off_qubits)
+
 from ._rdm_mapping_functions import (kronecker_delta,
                                      map_two_pdm_to_two_hole_dm,
                                      map_two_pdm_to_one_pdm,
diff --git a/src/openfermion/utils/_qubit_tapering_from_stabilizer.py b/src/openfermion/utils/_qubit_tapering_from_stabilizer.py
new file mode 100644
index 0000000..261c58e
--- /dev/null
+++ b/src/openfermion/utils/_qubit_tapering_from_stabilizer.py
@@ -0,0 +1,487 @@
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+
+"""Tools to reduce the number of terms and taper off qubits
+using stabilizer conditions. Based on ideas of arXiv:1701.08213. """
+
+import numpy
+from openfermion.ops import QubitOperator
+from openfermion.utils import count_qubits
+from openfermion.config import EQ_TOLERANCE
+
+
+class StabilizerError(Exception):
+    """Stabilizer error class."""
+
+    def __init__(self, message):
+        """
+        Throw custom errors connected to stabilizers.
+
+        Args:
+            message(str): custome error message string.
+        """
+        Exception.__init__(self, message)
+
+
+def check_commuting_stabilizers(stabilizer_list, msg, thres=EQ_TOLERANCE):
+    """
+    Auxiliary function checking that stabilizers commute.
+
+    If two stabilizers anti-commute their product
+    will have an imaginary coefficient.
+    This function checks the list of stabilizers (QubitOperator)
+    and raises and error if a complex number is found in
+    any of the coefficients.
+
+    Args:
+        stabilizer_list (list): List of stabilizers as QubitOperators.
+        msg (str): Message for the error.
+        thres: Tolerance value, set to OpenFermion tolerance by default.
+    """
+    for stab in stabilizer_list:
+        if abs(numpy.imag(list(stab.terms.values())[0])) >= thres:
+            raise StabilizerError(msg)
+
+
+def check_stabilizer_linearity(stabilizer_list, msg):
+    """
+    Auxiliary function checking that stabilizers are linearly independent.
+
+    If two stabilizer are linearly dependent the result
+    after some of their products will be the identity.
+    This function checks the list of stabilizers (QubitOperator)
+    and raises an error if the identity is found.
+
+    Args:
+        stabilizer_list (list): List of stabilizers (QubitOperator).
+        msg (str): Message for the error.
+    """
+    for stab in stabilizer_list:
+        if list(stab.terms.keys())[0] == ():
+            raise StabilizerError(msg)
+
+
+def fix_single_term(term, position, fixed_op, other_op, stabilizer):
+    """
+    Auxiliary function for term reductions.
+
+    Automatically multiplies a single term with a given stabilizer if
+    the Pauli operator on a given qubit is of one of two specified types.
+    This fixes a certain representation of a logical operator.
+
+    Args:
+        term (QubitOperator): Single term to fix.
+        position (int): Index of the qubit which is to be fixed.
+        fixed_op (str): Pauli operator, which
+                        will cause a multiplication by the
+                        stabilizer when encountered at the fixed
+                        position.
+        other_op (str): Alternative Pauli operator, which will also
+                        cause the multiplication by the stabilizer.
+        stabilizer (QubitOperator): Stabilizer that is multiplied
+                                    when necessary.
+
+    Returns:
+        term (QubitOperator): Updated term in a fiixed representation.
+    """
+    pauli_tuple = list(term.terms)[0]
+    if (position, fixed_op) in pauli_tuple or (position,
+                                               other_op) in pauli_tuple:
+        return term * stabilizer
+    else:
+        return term
+
+
+def _lookup_term(pauli_string, updated_terms_1, updated_terms_2):
+    """
+    Auxiliary function for reducing terms keeping length.
+
+    This function checks the length of the original Pauli strings,
+    compares it to a list of strings and returns the shortest operator
+    equivalent to the original.
+
+    Args:
+        pauli_string (tuple): Original Pauli string given in the same form
+                              as in the data structure of QubitOperator.
+        updated_terms_1 (list): List of Pauli strings (QubitOperator),
+                                which replace the original string if
+                                they are shorter and they are equivalent
+                                to each other.
+        updated_terms_2 (list): List of Pauli strings given in the data
+                                structure of QubitOperator, denoting which
+                                strings the entries of the first list are
+                                equivalent to.
+    Returns:
+        pauli_op (QubitOperator): Shortest Pauli string equivalent to the
+                                  original.
+
+    """
+    pauli_op = QubitOperator(pauli_string)
+    length = len(pauli_string)
+
+    for x in numpy.arange(len(updated_terms_1)):
+        if (pauli_string == updated_terms_2[x] and
+                (length > len(list(updated_terms_1[x].terms)[0]))):
+            pauli_op = updated_terms_1[x]
+            length = len(list(updated_terms_1[x].terms)[0])
+    return pauli_op
+
+
+def _reduce_terms(terms, stabilizer_list,
+                  manual_input, fixed_positions):
+    """
+    Perform the term reduction using stabilizer conditions.
+
+    Auxiliary function to reduce_number_of_terms.
+
+    Args:
+        terms (QubitOperator): Operator the number of terms is to be reduced.
+        stabilizer_list (list): List of the stabilizers as QubitOperators.
+        manual_input (Boolean): Option to pass the list of fixed qubits
+                                positions manually. Set to False by default.
+        fixed_positions (list): (optional) List of fixed qubit positions.
+                                Passing a list is only effective if
+                                manual_input is True.
+    Returns:
+        even_newer_terms (QubitOperator): Updated operator with reduced terms.
+        fixed_positions (list): Positions of qubits to be used for the
+                                term reduction.
+    Raises:
+        StabilizerError: Trivial stabilizer (identity).
+        StabilizerError: Stabilizer with complex coefficient.
+    """
+    # Initialize fixed_position as an empty list to avoid conflict with
+    # fixed_positions.
+    if manual_input is False:
+        fixed_positions = []
+
+    # We need the index of the stabilizer to connect it to the fixed qubit.
+    for i, x in enumerate(stabilizer_list):
+        selected_stab = list(stabilizer_list[0].terms)[0]
+
+        if manual_input is False:
+            # Find first position non-fixed position with non-trivial Pauli.
+            for qubit_pauli in selected_stab:
+                if qubit_pauli[0] not in fixed_positions:
+                    fixed_positions += [qubit_pauli[0]]
+                    fixed_op = qubit_pauli[1]
+                    break
+
+        else:
+            # Finds Pauli of the fixed qubit.
+            for qubit_pauli in selected_stab:
+                if qubit_pauli[0] == fixed_positions[i]:
+                    fixed_op = qubit_pauli[1]
+                    break
+
+        if fixed_op in ['X', 'Z']:
+            other_op = 'Y'
+        else:
+            other_op = 'X'
+
+        new_terms = QubitOperator()
+        for qubit_pauli in terms:
+            new_terms += fix_single_term(qubit_pauli, fixed_positions[i],
+                                         fixed_op,
+                                         other_op, stabilizer_list[0])
+        updated_stabilizers = []
+        for update_stab in stabilizer_list[1:]:
+            updated_stabilizers += [fix_single_term(update_stab,
+                                                    fixed_positions[i],
+                                                    fixed_op,
+                                                    other_op,
+                                                    stabilizer_list[0])]
+
+        # Update terms and stabilizer list.
+        terms = new_terms
+        stabilizer_list = updated_stabilizers
+
+        check_stabilizer_linearity(stabilizer_list,
+                                   msg='Linearly dependent stabilizers.')
+        check_commuting_stabilizers(stabilizer_list,
+                                    msg='Stabilizers anti-commute.')
+
+    return terms, fixed_positions
+
+
+def _reduce_terms_keep_length(terms, stabilizer_list,
+                              manual_input, fixed_positions):
+    """
+    Perform the term reduction using stabilizer conditions.
+
+    Auxiliary function to reduce_number_of_terms that returns the
+    Pauli strings with the same length as in the starting operator.
+
+    Args:
+        terms (QubitOperator): Operator from which terms are reduced.
+        stabilizer_list (list): List of the stabilizers as QubitOperators.
+        manual_input (Boolean): Option to pass the list of fixed qubits
+                                positions manually. Set to False by default.
+        fixed_positions (list): (optional) List of fixed qubit positions.
+                                Passing a list is only effective if
+                                manual_input is True.
+    Returns:
+        even_newer_terms (QubitOperator): Updated operator with reduced terms.
+        fixed_positions (list): Positions of qubits to be used for the
+                                term reduction.
+    Raises:
+        StabilizerError: Trivial stabilizer (identity).
+        StabilizerError: Stabilizer with complex coefficient.
+    """
+    term_list_duplicate = list(terms.terms)
+    term_list = [QubitOperator(x) for x in term_list_duplicate]
+
+    if manual_input is False:
+        fixed_positions = []
+
+    for i, x in enumerate(stabilizer_list):
+        selected_stab = list(stabilizer_list[0].terms)[0]
+
+        if manual_input is False:
+            # Finds qubit position and its Pauli.
+            for qubit_pauli in selected_stab:
+                if qubit_pauli[0] not in fixed_positions:
+                    fixed_positions += [qubit_pauli[0]]
+                    fixed_op = qubit_pauli[1]
+                    break
+        else:
+            # Finds Pauli of the fixed qubit.
+            for qubit_pauli in selected_stab:
+                if qubit_pauli[0] == fixed_positions[i]:
+                    fixed_op = qubit_pauli[1]
+                    break
+
+        if fixed_op in ['X', 'Z']:
+            other_op = 'Y'
+        else:
+            other_op = 'X'
+
+        new_list = []
+        updated_stabilizers = []
+        for y in term_list:
+            new_list += [fix_single_term(y, fixed_positions[i], fixed_op,
+                                         other_op, stabilizer_list[0])]
+        for update_stab in stabilizer_list[1:]:
+            updated_stabilizers += [fix_single_term(update_stab,
+                                                    fixed_positions[i],
+                                                    fixed_op, other_op,
+                                                    stabilizer_list[0])]
+        term_list = new_list
+        stabilizer_list = updated_stabilizers
+
+        check_stabilizer_linearity(stabilizer_list,
+                                   msg='Linearly dependent stabilizers.')
+        check_commuting_stabilizers(stabilizer_list,
+                                    msg='Stabilizers anti-commute.')
+
+    new_terms = QubitOperator()
+    for x, ent in enumerate(term_list):
+        new_terms += ent * terms.terms[term_list_duplicate[x]]
+    for x, ent in enumerate(term_list):
+        term_list_duplicate[x] = (
+            QubitOperator(term_list_duplicate[x]) /
+            list(ent.terms.items())[0][1])
+        term_list[x] = list(ent.terms)[0]
+
+    even_newer_terms = QubitOperator()
+    for pauli_string, coefficient in new_terms.terms.items():
+        even_newer_terms += coefficient * _lookup_term(
+            pauli_string,
+            term_list_duplicate,
+            term_list)
+
+    return even_newer_terms, fixed_positions
+
+
+def reduce_number_of_terms(operator, stabilizers,
+                           maintain_length=False,
+                           output_fixed_positions=False,
+                           manual_input=False,
+                           fixed_positions=None):
+    r"""
+    Reduce the number of Pauli strings of operator using stabilizers.
+
+    This function reduces the number of terms in a string by merging
+    terms that are identical by the multiplication of stabilizers.
+    The resulting Pauli strings maintain their length, unless specified
+    otherwise. In the latter case, a list of indices can be passed to
+    manually indicate the qubits to be fixed.
+
+    It is possible to reduce the number of terms in a Hamiltonian by
+    merging Pauli strings :math:`H_1, \, H_2` that are related by a
+    stabilizer :math:`S` such that  :math:`H_1 = H_2 \cdot S`. Given
+    a stabilizer generator :math:`\pm X \otimes p` this algorithm fixes the
+    first qubit, such that every Pauli string in the Hamiltonian acts with
+    either :math:`Z` or the identity on it. Where necessary, this is achieved
+    by multiplications with :math:`\pm X \otimes p`: a string
+    :math:`Y \otimes h`, for instance,  is turned into
+    :math:`Z \otimes (\mp ih\cdot p)`. Qubits on which a generator acts as
+    :math:`Y` (:math:`Z`) are constrained to be acted on by the Hamiltonian as
+    :math:`Z` (:math:`X`) or the identity. Fixing a different qubit for every
+    stabilizer generator eliminates all redundant strings. The fixed
+    representations are in the end re-expressed as the shortest of the
+    original strings, :math:`H_1` or :math:`H_2`.
+
+
+    Args:
+        operator (QubitOperator): Operator of which the number of terms
+                                  will be reduced.
+        stabilizers (QubitOperator): Stabilizer generators used for the
+                                     reduction. Can also be passed as a list
+                                     of QubitOperator.
+        maintain_length (Boolean): Option deciding whether the fixed Pauli
+                                   strings are re-expressed in their original
+                                   form. Set to False by default.
+        output_fixed_positions (Boolean): Option deciding whether to return
+                                          the list of fixed qubit positions.
+                                          Set to False by default.
+        manual_input (Boolean): Option to pass the list of fixed qubits
+                                positions manually. Set to False by default.
+        fixed_positions (list): (optional) List of fixed qubit positions.
+                                Passing a list is only effective if
+                                manual_input is True.
+    Returns:
+        reduced_operator (QubitOperator): Operator with reduced number of
+                                          terms.
+
+        fixed_positions (list): (optional) Fixed qubits.
+
+    Raises:
+        TypeError: Input terms must be QubitOperator.
+        TypeError: Input stabilizers must be QubitOperator or list.
+        StabilizerError: Trivial stabilizer (identity).
+        StabilizerError: Stabilizer with complex coefficient.
+        TypeError: List of fixed qubits required if manual input True.
+        StabilizerError: The number of stabilizers must be equal to the number
+                         of qubits manually fixed.
+        StabilizerError: All qubit positions must be different.
+    """
+    if not isinstance(operator, QubitOperator):
+        raise TypeError('Input terms must be QubitOperator.')
+    if not isinstance(stabilizers, (QubitOperator, list,
+                                    tuple, numpy.ndarray)):
+        raise TypeError('Input stabilizers must be QubitOperator or list.')
+
+    stabilizer_list = list(stabilizers)
+
+    check_stabilizer_linearity(stabilizer_list,
+                               msg='Trivial stabilizer (identity).')
+    check_commuting_stabilizers(stabilizer_list,
+                                msg='Stabilizer with complex coefficient.')
+
+    if manual_input:
+        # Convert fixed_position into a list to allow any type of
+        # array_like data structure.
+        if fixed_positions is None:
+            raise TypeError('List of qubit positions required.')
+        fixed_positions = list(fixed_positions)
+        if len(fixed_positions) != len(stabilizer_list):
+            raise StabilizerError('The number of stabilizers must be equal '
+                                  + 'to the number of qubits manually fixed.')
+        if len(set(fixed_positions)) != len(stabilizer_list):
+            raise StabilizerError('All qubit positions must be different.')
+
+    if maintain_length:
+        (reduced_operator,
+         fixed_positions) = _reduce_terms_keep_length(operator,
+                                                      stabilizer_list,
+                                                      manual_input,
+                                                      fixed_positions)
+    else:
+        (reduced_operator,
+         fixed_positions) = _reduce_terms(operator,
+                                          stabilizer_list,
+                                          manual_input,
+                                          fixed_positions)
+
+    if output_fixed_positions:
+        return reduced_operator, fixed_positions
+    else:
+        return reduced_operator
+
+
+def taper_off_qubits(operator, stabilizers, manual_input=False,
+                     fixed_positions=None, output_tapered_positions=False):
+    r"""
+    Remove qubits from given operator.
+
+    Qubits are removed by eliminating an equivalent number of
+    stabilizer conditions. Which qubits that are can either be determined
+    automatically or their positions can be set manually.
+
+    Qubits can be disregarded from the Hamiltonian when the effect of all its
+    terms on them is rendered trivial. This algorithm employs a stabilizers
+    like :math:`\pm X \otimes p` to fix the action of every Pauli
+    string on the first qubit to :math:`Z` or the identity. A string
+    :math:`X \otimes h` would for instance be multiplied with the stabilizer
+    to obtain :math:`1 \otimes (\pm h\cdot p)` while a string
+    :math:`Z \otimes h^\prime` would pass without correction. The first
+    qubit can subsequently be removed as it must be in the computational basis
+    in Hamiltonian eigenstates.
+    For stabilizers acting as :math:`Y` (:math:`Z`) on selected qubits,
+    the algorithm would fix the action of every Hamiltonian string to
+    :math:`Z` (:math:`X`). Updating also the list of remaining stabilizer
+    generators, the algorithm is run iteratively.
+
+    Args:
+        operator (QubitOperator): Operator of which qubits will be removed.
+        stabilizers (QubitOperator): Stabilizer generators for the tapering.
+                                     Can also be passed as a list of
+                                     QubitOperator.
+        manual_input (Boolean): Option to pass the list of fixed qubits
+                                positions manually. Set to False by default.
+        fixed_positions (list): (optional) List of fixed qubit positions.
+                                Passing a list is only effective if
+                                manual_input is True.
+        output_tapered_positions (Boolean): Option to output the positions of
+                                            qubits that have been removed.
+    Returns:
+        skimmed_operator (QubitOperator): Operator with fewer qubits.
+        removed_positions (list): (optional) List of removed qubit positions.
+                                  For the qubits to be gone in the qubit count,
+                                  the remaining qubits have been moved up to
+                                  those indices.
+    """
+    n_qbits = count_qubits(operator)
+    (ham_to_update,
+     qbts_to_rm) = reduce_number_of_terms(operator,
+                                          stabilizers,
+                                          maintain_length=False,
+                                          manual_input=manual_input,
+                                          fixed_positions=fixed_positions,
+                                          output_fixed_positions=True)
+
+    # Gets a list of the order of the qubits after tapering
+    qbit_order = list(numpy.arange(n_qbits - len(qbts_to_rm), dtype=int))
+    # Save the original list before it gets ordered
+    removed_positions = qbts_to_rm
+    qbts_to_rm.sort()
+    for x in qbts_to_rm:
+        qbit_order.insert(x, 'remove')
+
+    # Remove the qubits
+    skimmed_operator = QubitOperator()
+    for term, coef in ham_to_update.terms.items():
+        if term == ():
+            skimmed_operator += QubitOperator('', coef)
+            continue
+        tap_tpls = []
+        for p in term:
+            if qbit_order[p[0]] != 'remove':
+                tap_tpls.append((qbit_order[p[0]].item(), p[1]))
+
+        skimmed_operator += QubitOperator(tuple(tap_tpls), coef)
+
+    if output_tapered_positions:
+        return skimmed_operator, removed_positions
+    else:
+        return skimmed_operator


Test patch:
diff --git a/src/openfermion/utils/_qubit_tapering_from_stabilizer_test.py b/src/openfermion/utils/_qubit_tapering_from_stabilizer_test.py
new file mode 100644
index 0000000..90c0da4
--- /dev/null
+++ b/src/openfermion/utils/_qubit_tapering_from_stabilizer_test.py
@@ -0,0 +1,337 @@
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+
+"""Test for qubit_tapering_from_stabilizer model."""
+
+import unittest
+
+import numpy
+
+from openfermion.hamiltonians import MolecularData
+from openfermion.ops import QubitOperator
+from openfermion.transforms import jordan_wigner, get_fermion_operator
+from openfermion.utils import eigenspectrum
+
+from openfermion.utils import reduce_number_of_terms, taper_off_qubits
+from openfermion.utils._qubit_tapering_from_stabilizer import (
+    StabilizerError, check_commuting_stabilizers, check_stabilizer_linearity,
+    fix_single_term, _reduce_terms, _reduce_terms_keep_length,
+    _lookup_term)
+
+
+def lih_hamiltonian():
+    """
+    Generate test Hamiltonian from LiH.
+
+    Args:
+        None
+
+    Return:
+
+        hamiltonian: FermionicOperator
+
+        spectrum: List of energies.
+    """
+    geometry = [('Li', (0., 0., 0.)), ('H', (0., 0., 1.45))]
+    active_space_start = 1
+    active_space_stop = 3
+    molecule = MolecularData(geometry, 'sto-3g', 1,
+                             description="1.45")
+    molecule.load()
+
+    molecular_hamiltonian = molecule.get_molecular_hamiltonian(
+        occupied_indices=range(active_space_start),
+        active_indices=range(active_space_start, active_space_stop))
+
+    hamiltonian = get_fermion_operator(molecular_hamiltonian)
+    spectrum = eigenspectrum(hamiltonian)
+
+    return hamiltonian, spectrum
+
+
+class TaperingTest(unittest.TestCase):
+    """TaperingTest class."""
+
+    def test_function_errors(self):
+        """Test error of main function."""
+        hamiltonian, spectrum = lih_hamiltonian()
+        qubit_hamiltonian = jordan_wigner(hamiltonian)
+        stab1 = QubitOperator('Z0 Z2', -1.0)
+        stab2 = QubitOperator('Z1 Z3', -1.0)
+        with self.assertRaises(TypeError):
+            reduce_number_of_terms(operator=1,
+                                   stabilizers=stab1 + stab2)
+        with self.assertRaises(TypeError):
+            reduce_number_of_terms(operator=qubit_hamiltonian,
+                                   stabilizers=1)
+        with self.assertRaises(TypeError):
+            reduce_number_of_terms(operator=qubit_hamiltonian,
+                                   stabilizers=stab1 + stab2,
+                                   manual_input=True,
+                                   fixed_positions=None)
+        with self.assertRaises(StabilizerError):
+            reduce_number_of_terms(operator=qubit_hamiltonian,
+                                   stabilizers=stab1 + stab2,
+                                   manual_input=True,
+                                   fixed_positions=[1])
+        with self.assertRaises(StabilizerError):
+            reduce_number_of_terms(operator=qubit_hamiltonian,
+                                   stabilizers=stab1 + stab2,
+                                   manual_input=True,
+                                   fixed_positions=[1, 1])
+        with self.assertRaises(StabilizerError):
+            # Check Identity as stabilizer error.
+            reduce_number_of_terms(operator=qubit_hamiltonian,
+                                   stabilizers=(stab1 +
+                                                QubitOperator(' ', 1.0)))
+        with self.assertRaises(StabilizerError):
+            # Check complex coefficient stabilizer error.
+            reduce_number_of_terms(operator=qubit_hamiltonian,
+                                   stabilizers=(stab1 +
+                                                QubitOperator('Z0', 1.0j)))
+        with self.assertRaises(StabilizerError):
+            # Check linearly-dependent stabilizer error.
+            reduce_number_of_terms(operator=qubit_hamiltonian,
+                                   stabilizers=(stab1 +
+                                                QubitOperator('Z0 Z1 Z2 Z3',
+                                                              1.0) +
+                                                stab2))
+        with self.assertRaises(StabilizerError):
+            # Check anti-commuting stabilizer error.
+            reduce_number_of_terms(operator=qubit_hamiltonian,
+                                   stabilizers=(QubitOperator('X0', 1.0) +
+                                                QubitOperator('Y0', 1.0)))
+        with self.assertRaises(StabilizerError):
+            # Check linearly-dependent stabilizer error.
+            _reduce_terms(terms=qubit_hamiltonian,
+                          stabilizer_list=list(stab1 +
+                                               QubitOperator('Z0 Z1 Z2 Z3',
+                                                             1.0) +
+                                               stab2),
+                          manual_input=False,
+                          fixed_positions=[])
+        with self.assertRaises(StabilizerError):
+            # Check complex coefficient stabilizer error.
+            _reduce_terms(terms=qubit_hamiltonian,
+                          stabilizer_list=list(stab1 +
+                                               QubitOperator('Z0', 1.0j)),
+                          manual_input=False,
+                          fixed_positions=[])
+        with self.assertRaises(StabilizerError):
+            # Check linearly-dependent stabilizer error.
+            par_qop = QubitOperator('Z0 Z1 Z2 Z3', 1.0)
+            _reduce_terms_keep_length(terms=qubit_hamiltonian,
+                                      stabilizer_list=[stab1,
+                                                       par_qop,
+                                                       stab2],
+                                      manual_input=False,
+                                      fixed_positions=[])
+        with self.assertRaises(StabilizerError):
+            # Check complex coefficient stabilizer error.
+            aux_qop = QubitOperator('Z0',
+                                    1.0j)
+            _reduce_terms_keep_length(terms=qubit_hamiltonian,
+                                      stabilizer_list=[stab1,
+                                                       aux_qop],
+                                      manual_input=False,
+                                      fixed_positions=[])
+        with self.assertRaises(StabilizerError):
+            # Test check_commuting_stabilizer function
+            # Requires a list of QubitOperators one of which
+            # has an imaginary term.
+            check_commuting_stabilizers(stabilizer_list=[QubitOperator('Z0 Z1',
+                                                                       1.0),
+                                                         QubitOperator('X0',
+                                                                       1j)],
+                                        msg='This test fails.')
+        with self.assertRaises(StabilizerError):
+            # Test check_stabilizer_linearity function.
+            # Requires a list of QUbitOperators one of which is
+            # the identity.
+            check_stabilizer_linearity([QubitOperator('Z0 Z1', 1.0),
+                                        QubitOperator(' ', 1.0)],
+                                       msg='This test fails.')
+
+    def test_fix_single_term(self):
+        """Test fix_single_term function."""
+        stab2 = QubitOperator('Z1 Z3', -1.0)
+        test_term = QubitOperator('Z1 Z2')
+
+        fix1 = fix_single_term(test_term, 1, 'Z', 'X', stab2)
+        fix2 = fix_single_term(test_term, 0, 'X', 'X', stab2)
+
+        self.assertTrue(fix1 == (test_term * stab2))
+        self.assertTrue(fix2 == test_term)
+
+    def test_lookup_term(self):
+        """Test for the auxiliar function _lookup_term."""
+        # Dummy test where the initial Pauli string is larger than the
+        # updated one.
+        start_op = list(QubitOperator('Z0 Z1 Z2 Z3').terms.keys())[0]
+        updateop1 = QubitOperator('Z0 Z2', -1.0)
+        updateop2 = list(QubitOperator('Z0 Z1 Z2 Z3').terms.keys())
+
+        qop = _lookup_term(start_op, [updateop1], updateop2)
+        final_op = list(qop.terms.keys())[0]
+
+        self.assertLess(len(final_op), len(start_op))
+
+    def test_reduce_terms(self):
+        """Test reduce_terms function using LiH Hamiltonian."""
+        hamiltonian, spectrum = lih_hamiltonian()
+        qubit_hamiltonian = jordan_wigner(hamiltonian)
+        stab1 = QubitOperator('Z0 Z2', -1.0)
+        stab2 = QubitOperator('Z1 Z3', -1.0)
+
+        red_eigenspectrum = eigenspectrum(
+            reduce_number_of_terms(qubit_hamiltonian,
+                                   stab1 + stab2))
+
+        self.assertAlmostEqual(spectrum[0], red_eigenspectrum[0])
+
+    def test_reduce_terms_manual_input(self):
+        """Test reduce_terms function using LiH Hamiltonian."""
+        hamiltonian, spectrum = lih_hamiltonian()
+        qubit_hamiltonian = jordan_wigner(hamiltonian)
+        stab1 = QubitOperator('Z0 Z2', -1.0)
+        stab2 = QubitOperator('Z1 Z3', -1.0)
+
+        red_eigenspectrum = eigenspectrum(
+            reduce_number_of_terms(qubit_hamiltonian,
+                                   [stab1, stab2],
+                                   manual_input=True,
+                                   fixed_positions=[0, 1]))
+
+        self.assertAlmostEqual(spectrum[0], red_eigenspectrum[0])
+
+    def test_reduce_terms_maintain_length(self):
+        """Test reduce_terms function using LiH Hamiltonian."""
+        hamiltonian, spectrum = lih_hamiltonian()
+        qubit_hamiltonian = jordan_wigner(hamiltonian)
+        stab1 = QubitOperator('Z0 Z2', -1.0)
+        stab2 = QubitOperator('Z1 Z3', -1.0)
+
+        red_eigenspectrum = eigenspectrum(
+            reduce_number_of_terms(qubit_hamiltonian,
+                                   stab1 + stab2,
+                                   maintain_length=True))
+
+        self.assertAlmostEqual(spectrum[0], red_eigenspectrum[0])
+
+    def test_reduce_terms_auxiliar_functions(self):
+        """Test reduce_terms function using LiH Hamiltonian."""
+        hamiltonian, spectrum = lih_hamiltonian()
+        qubit_ham = jordan_wigner(hamiltonian)
+        stab1 = QubitOperator('Z0 Z2', -1.0)
+        stab2 = QubitOperator('Z1 Z3', -1.0)
+
+        red_ham1, fixed_pos1 = _reduce_terms(terms=qubit_ham,
+                                             stabilizer_list=[stab1, stab2],
+                                             manual_input=False,
+                                             fixed_positions=[])
+        red_ham2, fixed_pos2 = _reduce_terms_keep_length(terms=qubit_ham,
+                                                         stabilizer_list=[
+                                                             stab1, stab2],
+                                                         manual_input=False,
+                                                         fixed_positions=[])
+        red_eigspct1 = eigenspectrum(red_ham1)
+        red_eigspct2 = eigenspectrum(red_ham2)
+
+        self.assertAlmostEqual(spectrum[0], red_eigspct1[0])
+        self.assertAlmostEqual(spectrum[0], red_eigspct2[0])
+
+    def test_reduce_terms_auxiliar_functions_manual_input(self):
+        """Test reduce_terms function using LiH Hamiltonian."""
+        hamiltonian, spectrum = lih_hamiltonian()
+        qubit_ham = jordan_wigner(hamiltonian)
+        stab1 = QubitOperator('Z0 Z2', -1.0)
+        stab2 = QubitOperator('Z1 Z3', -1.0)
+
+        red_ham1, fixed_pos1 = _reduce_terms(terms=qubit_ham,
+                                             stabilizer_list=[stab1, stab2],
+                                             manual_input=True,
+                                             fixed_positions=[0, 1])
+        red_ham2, fixed_pos2 = _reduce_terms_keep_length(terms=qubit_ham,
+                                                         stabilizer_list=[
+                                                             stab1, stab2],
+                                                         manual_input=True,
+                                                         fixed_positions=[0,
+                                                                          1])
+        red_eigspct1 = eigenspectrum(red_ham1)
+        red_eigspct2 = eigenspectrum(red_ham2)
+
+        self.assertAlmostEqual(spectrum[0], red_eigspct1[0])
+        self.assertAlmostEqual(spectrum[0], red_eigspct2[0])
+
+    def test_tapering_qubits_manual_input_false(self):
+        """Test taper_off_qubits function using LiH Hamiltonian."""
+        hamiltonian, spectrum = lih_hamiltonian()
+        qubit_hamiltonian = jordan_wigner(hamiltonian)
+        stab1 = QubitOperator('Z0 Z2', -1.0)
+        stab2 = QubitOperator('Z1 Z3', -1.0)
+
+        tapered_hamiltonian = taper_off_qubits(operator=qubit_hamiltonian,
+                                               stabilizers=[stab1, stab2],
+                                               manual_input=False,
+                                               fixed_positions=[0, 3])
+        tapered_spectrum = eigenspectrum(tapered_hamiltonian)
+
+        self.assertAlmostEqual(spectrum[0], tapered_spectrum[0])
+
+    def test_tapering_qubits_manual_input(self):
+        """
+        Test taper_off_qubits function using LiH Hamiltonian.
+
+        Checks different qubits inputs to remove manually.
+
+        Test the lowest eigenvalue against the full Hamiltonian,
+        and the full spectrum between them.
+        """
+        hamiltonian, spectrum = lih_hamiltonian()
+        qubit_hamiltonian = jordan_wigner(hamiltonian)
+        stab1 = QubitOperator('Z0 Z2', -1.0)
+        stab2 = QubitOperator('Z1 Z3', -1.0)
+
+        tapered_ham_0_3 = taper_off_qubits(qubit_hamiltonian,
+                                           [stab1, stab2],
+                                           manual_input=True,
+                                           fixed_positions=[0, 3])
+        tapered_ham_2_1 = taper_off_qubits(qubit_hamiltonian,
+                                           [stab1, stab2],
+                                           manual_input=True,
+                                           fixed_positions=[2, 1])
+
+        tapered_spectrum_0_3 = eigenspectrum(tapered_ham_0_3)
+        tapered_spectrum_2_1 = eigenspectrum(tapered_ham_2_1)
+
+        self.assertAlmostEqual(spectrum[0], tapered_spectrum_0_3[0])
+        self.assertAlmostEqual(spectrum[0], tapered_spectrum_2_1[0])
+        self.assertTrue(numpy.allclose(tapered_spectrum_0_3,
+                                       tapered_spectrum_2_1))
+
+    def test_tapering_qubits_remove_positions(self):
+        """Test taper_off_qubits function using LiH Hamiltonian."""
+        hamiltonian, spectrum = lih_hamiltonian()
+        qubit_hamiltonian = jordan_wigner(hamiltonian)
+        stab1 = QubitOperator('Z0 Z2', -1.0)
+        stab2 = QubitOperator('Z1 Z3', -1.0)
+
+        (tapered_hamiltonian,
+         positions) = taper_off_qubits(operator=qubit_hamiltonian,
+                                       stabilizers=[stab1, stab2],
+                                       manual_input=True,
+                                       fixed_positions=[0, 3],
+                                       output_tapered_positions=True)
+        tapered_spectrum = eigenspectrum(tapered_hamiltonian)
+
+        self.assertAlmostEqual(spectrum[0], tapered_spectrum[0])
+        self.assertEqual(positions, [0, 3])