Repository: quantumlib/OpenFermion
Pull request #559 (2020-02-26 20:44:55)

Remove trailing whitepaces from .py files
Part of https://github.com/quantumlib/OpenFermion/issues/416

vtomole: @kevinsung I haven't added or removed code. All I did was format. Is this one of those false positives you were talking about?

Patch:
diff --git a/src/openfermion/hamiltonians/_general_hubbard.py b/src/openfermion/hamiltonians/_general_hubbard.py
index 0490fee..2da88b4 100644
--- a/src/openfermion/hamiltonians/_general_hubbard.py
+++ b/src/openfermion/hamiltonians/_general_hubbard.py
@@ -35,14 +35,15 @@ def number_operator(i, coefficient=1., particle_hole_symmetry=False):
 
 
 def interaction_operator(i, j, coefficient=1., particle_hole_symmetry=False):
-    return (number_operator(i, coefficient, 
-                            particle_hole_symmetry=particle_hole_symmetry) * 
+    return (number_operator(
+        i, coefficient, particle_hole_symmetry=particle_hole_symmetry) *
             number_operator(j, particle_hole_symmetry=particle_hole_symmetry))
 
 
 def tunneling_operator(i, j, coefficient=1.):
-    return (FermionOperator(((i, 1), (j, 0)), coefficient) + 
-            FermionOperator(((j, 1), (i, 0)), coefficient.conjugate()))
+    return (FermionOperator(((i, 1), (j, 0)), coefficient) + FermionOperator(
+        ((j, 1), (i, 0)), coefficient.conjugate()))
+
 
 def number_difference_operator(i, j, coefficient=1.):
     return number_operator(i, coefficient) - number_operator(j, coefficient)
@@ -195,20 +196,20 @@ class FermiHubbardModel:
 
         .. math::
 
-            t \sum_{(i, j) \in E^{(\mathrm{edge type})}} 
+            t \sum_{(i, j) \in E^{(\mathrm{edge type})}}
             \sum_{\sigma}
-            \left(a_{i, a, \sigma}^{\dagger} a_{j, b, \sigma} 
+            \left(a_{i, a, \sigma}^{\dagger} a_{j, b, \sigma}
             + a_{j, b, \sigma}^{\dagger} a_{i, a, \sigma}\right)
 
-        and in the spinless model to 
+        and in the spinless model to
 
         .. math::
 
-            -t \sum_{(i, j) \in E^{(\mathrm{edge type})}} 
-            \left(a_{i, a}^{\dagger} a_{j, b} 
+            -t \sum_{(i, j) \in E^{(\mathrm{edge type})}}
+            \left(a_{i, a}^{\dagger} a_{j, b}
             + a_{j, b}^{\dagger} a_{i, a}\right),
 
-        where 
+        where
 
             - :math:`(a, b)` is the pair of degrees of freedom given by ``dofs``;
             - :math:`E^{(\mathrm{edge type})}` is the set of ordered pairs of
@@ -228,7 +229,7 @@ class FermiHubbardModel:
             U \sum_{(i, j) \in E^{(\mathrm{edge type})}} \sum_{(\sigma, \sigma')}
             n_{i, a, \sigma} n_{j, b, \sigma'}
 
-        where 
+        where
 
             - :math:`(a, b)` is the pair of degrees of freedom given by ``dofs``;
             - :math:`E^{(\mathrm{edge type})}` is the set of ordered pairs of
@@ -236,7 +237,7 @@ class FermiHubbardModel:
             - :math:`U` is the ``coefficient``; and
             - :math:`(\sigma, \sigma')` runs over
                 - all four possible pairs of spins if `spin_pairs == SpinPairs.ALL`,
-                - :math:`\{(\uparrow, \downarrow), (\downarrow, \uparrow)\}` if `spin_pairs == SpinPairs.DIFF`, and 
+                - :math:`\{(\uparrow, \downarrow), (\downarrow, \uparrow)\}` if `spin_pairs == SpinPairs.DIFF`, and
                 - :math:`\{(\uparrow, \uparrow), (\downarrow, \downarrow)\}' if 'spin_pairs == SpinPairs.SAME`.
 
         Each potential parameter is a tuple ``(dof, coefficient)``. For example, in the spinful model, it corresponds to the terms
@@ -244,7 +245,7 @@ class FermiHubbardModel:
         .. math::
             -\mu \sum_{i} \sum_{\sigma} n_{i, a, \sigma},
 
-        where 
+        where
 
             - :math:`i` runs over the sites of the lattice;
             - :math:`a` is the degree of freedom ``dof``; and
@@ -273,7 +274,7 @@ class FermiHubbardModel:
             parameter = TunnelingParameter(*parameter)
             self.lattice.validate_edge_type(parameter.edge_type)
             self.lattice.validate_dofs(parameter.dofs, 2)
-            if ((parameter.edge_type in self.lattice.onsite_edge_types) and 
+            if ((parameter.edge_type in self.lattice.onsite_edge_types) and
                 (len(set(parameter.dofs)) == 1)):
                 raise ValueError('Invalid onsite tunneling parameter between '
                                  'same dof {}.'.format(parameter.dofs))
@@ -288,13 +289,13 @@ class FermiHubbardModel:
         for parameter in parameters:
             if len(parameter) not in (3, 4):
                 raise ValueError('len(parameter) not in (3, 4)')
-            spin_pairs = (SpinPairs.ALL if len(parameter) < 4 
-                                        else parameter[-1])
+            spin_pairs = (SpinPairs.ALL
+                          if len(parameter) < 4 else parameter[-1])
             parameter = InteractionParameter(*parameter[:3],
                                              spin_pairs=spin_pairs)
             self.lattice.validate_edge_type(parameter.edge_type)
             self.lattice.validate_dofs(parameter.dofs, 2)
-            if ((len(set(parameter.dofs)) == 1) and 
+            if ((len(set(parameter.dofs)) == 1) and
                 (parameter.edge_type in self.lattice.onsite_edge_types) and
                 (parameter.spin_pairs == SpinPairs.SAME)):
                 raise ValueError('Parameter {} specifies '.format(parameter) +
@@ -351,8 +352,10 @@ class FermiHubbardModel:
                 for spin_index in self.lattice.spin_indices:
                     i = self.lattice.to_spin_orbital_index(
                             site_index, param.dof, spin_index)
-                    terms += number_operator(i, -param.coefficient, 
-                            particle_hole_symmetry=self.particle_hole_symmetry)
+                    terms += number_operator(
+                        i,
+                        -param.coefficient,
+                        particle_hole_symmetry=self.particle_hole_symmetry)
         return terms
 
 
@@ -369,7 +372,5 @@ class FermiHubbardModel:
 
 
     def hamiltonian(self):
-        return (self.tunneling_terms() + 
-                self.interaction_terms() + 
-                self.potential_terms() +
-                self.field_terms())
+        return (self.tunneling_terms() + self.interaction_terms() +
+                self.potential_terms() + self.field_terms())
diff --git a/src/openfermion/hamiltonians/_molecular_data.py b/src/openfermion/hamiltonians/_molecular_data.py
index 7d32199..965026f 100644
--- a/src/openfermion/hamiltonians/_molecular_data.py
+++ b/src/openfermion/hamiltonians/_molecular_data.py
@@ -972,7 +972,7 @@ def load_molecular_hamiltonian(
     else:
         n_core_orbitals = (molecule.n_electrons - n_active_electrons) // 2
         occupied_indices = list(range(n_core_orbitals))
-        
+
     if n_active_orbitals is None:
         active_indices = None
     else:
diff --git a/src/openfermion/measurements/_qubit_partitioning.py b/src/openfermion/measurements/_qubit_partitioning.py
index a23a070..79c3b5d 100644
--- a/src/openfermion/measurements/_qubit_partitioning.py
+++ b/src/openfermion/measurements/_qubit_partitioning.py
@@ -67,7 +67,7 @@ def binary_partition_iterator(qubit_list, num_iterations=None):
         # which we delete.
         if qubit_list[-1] is None:
             del qubit_list[-1]
-        
+
 
 def partition_iterator(qubit_list, partition_size, num_iterations=None):
     """Generator for a list of k-partitions of N qubits such that
diff --git a/src/openfermion/ops/_majorana_operator.py b/src/openfermion/ops/_majorana_operator.py
index 8652c17..5fe7c52 100644
--- a/src/openfermion/ops/_majorana_operator.py
+++ b/src/openfermion/ops/_majorana_operator.py
@@ -96,7 +96,7 @@ class MajoranaOperator:
         The input to this method is a real orthogonal matrix :math:`O`.
         It returns a new MajoranaOperator which is equivalent to the old one
         but rewritten in terms of a new basis of Majorana operators.
-        Let the original Majorana operators be denoted by 
+        Let the original Majorana operators be denoted by
         :math:`\gamma_i` and the new operators be denoted by
         :math:`\tilde{\gamma_i}`. Then they are related by the equation
 
diff --git a/src/openfermion/ops/_polynomial_tensor.py b/src/openfermion/ops/_polynomial_tensor.py
index ea926b2..82c8137 100644
--- a/src/openfermion/ops/_polynomial_tensor.py
+++ b/src/openfermion/ops/_polynomial_tensor.py
@@ -362,7 +362,7 @@ class PolynomialTensor(object):
             selection (Union[int, Iterable[int]): If int, keeps terms with at
                 most (exactly, if exact is True) that many unique indices. If
                 iterable, keeps only terms containing (all of, if exact is
-                True) the specified indices. 
+                True) the specified indices.
             exact (bool): Whether or not the selection is strict.
         """
         comparator = (operator.eq if exact else operator.le)
diff --git a/src/openfermion/transforms/_binary_code_transform.py b/src/openfermion/transforms/_binary_code_transform.py
index 3dcdfb3..f313bee 100644
--- a/src/openfermion/transforms/_binary_code_transform.py
+++ b/src/openfermion/transforms/_binary_code_transform.py
@@ -71,7 +71,7 @@ def dissolve(term):
 
 def make_parity_list(code):
     """Create the parity list from the decoder of the input code.
-    The output parity list has a similar structure as code.decoder. 
+    The output parity list has a similar structure as code.decoder.
 
     Args:
         code (BinaryCode): the code to extract the parity list from.
@@ -102,7 +102,7 @@ def binary_code_transform(hamiltonian, code):
     property N>n, when the Fermion basis is smaller than the fermionic Fock
     space. The binary_code_transform function can transform Fermion operators
     to qubit operators for custom- and qubit-saving mappings.
-    
+
     Note:
         Logic multi-qubit operators are decomposed into Pauli-strings (e.g.
         CPhase(1,2) = 0.5 * (1 + Z1 + Z2 - Z1 Z2 ) ), which might increase
diff --git a/src/openfermion/transforms/_binary_codes.py b/src/openfermion/transforms/_binary_codes.py
index cc7726f..3d7012b 100644
--- a/src/openfermion/transforms/_binary_codes.py
+++ b/src/openfermion/transforms/_binary_codes.py
@@ -139,7 +139,7 @@ def checksum_code(n_modes, odd):
     is defined such that it yields the total occupation number for a given basis
     state. A Checksum code with odd weight will encode all states with odd
     occupation number. This code saves one qubit: n_qubits = n_modes - 1.
-        
+
     Args:
         n_modes (int): number of modes
         odd (int or bool): 1 (True) or 0 (False), if odd,
@@ -153,7 +153,7 @@ def checksum_code(n_modes, odd):
 
 def jordan_wigner_code(n_modes):
     """ The Jordan-Wigner transform as binary code.
-        
+
     Args:
         n_modes (int): number of modes
 
@@ -166,7 +166,7 @@ def jordan_wigner_code(n_modes):
 def bravyi_kitaev_code(n_modes):
     """ The Bravyi-Kitaev transform as binary code. The implementation
     follows arXiv:1208.5986.
-        
+
     Args:
         n_modes (int): number of modes
 
@@ -180,7 +180,7 @@ def parity_code(n_modes):
     """ The parity transform (arXiv:1208.5986) as binary code. This code is
     very similar to the Jordan-Wigner transform, but with long update strings
     instead of parity strings. It does not save qubits: n_qubits = n_modes.
-        
+
     Args:
         n_modes (int): number of modes
 
@@ -200,7 +200,7 @@ def weight_one_binary_addressing_code(exponent):
     of two. It encodes all possible vectors with Hamming weight 1, which
     corresponds to all states with total occupation one. The amount of
     qubits saved here is maximal: for a given argument 'exponent', we find
-    n_modes = 2 ^ exponent, n_qubits = exponent. 
+    n_modes = 2 ^ exponent, n_qubits = exponent.
 
     Note:
         This code is highly non-linear and might produce a lot of terms.
@@ -258,16 +258,16 @@ def weight_two_segment_code():
 def interleaved_code(modes):
     """ Linear code that reorders orbitals from even-odd to up-then-down.
     In up-then-down convention, one can append two instances of the same
-    code 'c' in order to have two symmetric subcodes that are symmetric for 
+    code 'c' in order to have two symmetric subcodes that are symmetric for
     spin-up and -down modes: ' c + c '.
-    In even-odd, one can concatenate with the interleaved_code 
+    In even-odd, one can concatenate with the interleaved_code
     to have the same result:' interleaved_code * (c + c)'.
     This code changes the order of modes from (0, 1 , 2, ... , modes-1 )
-    to (0, modes/2, 1 modes/2+1, ... , modes-1, modes/2 - 1).  
-    n_qubits = n_modes. 
-    
-    Args: modes (int): number of modes, must be even 
-    
+    to (0, modes/2, 1 modes/2+1, ... , modes-1, modes/2 - 1).
+    n_qubits = n_modes.
+
+    Args: modes (int): number of modes, must be even
+
     Returns (BinaryCode): code that interleaves orbitals
     """
     if modes % 2 == 1:
diff --git a/src/openfermion/transforms/_bksf.py b/src/openfermion/transforms/_bksf.py
index f47e150..7adbc81 100644
--- a/src/openfermion/transforms/_bksf.py
+++ b/src/openfermion/transforms/_bksf.py
@@ -68,7 +68,7 @@ def bravyi_kitaev_fast_interaction_op(iop):
     For the sake of brevity the reader is encouraged to look up the
     expressions of other terms from the code below. The variables for edge
     operators are chosen according to the nomenclature defined above
-    (B_i and A_ij). A detailed description of these operators and the terms 
+    (B_i and A_ij). A detailed description of these operators and the terms
     of the electronic Hamiltonian are provided in (arXiv 1712.00446).
 
     Args:
@@ -116,7 +116,7 @@ def bravyi_kitaev_fast_interaction_op(iop):
                             qubit_operator += transformed_term
                             continue
                         elif p != r and q < p:
-                                continue
+                            continue
 
 
                     # Handle the two-body terms.
@@ -167,7 +167,7 @@ def bravyi_kitaev_fast_edge_matrix(iop, n_qubits=None):
                             if min(r, s) < min(p, q):
                                 continue
                         elif p != r and q < p:
-                                continue
+                            continue
 
                     # Handle case of four unique indices.
                     if len(set([p, q, r, s])) == 4:
@@ -365,7 +365,7 @@ def edge_operator_aij(edge_matrix_indices, i, j):
     for edge_index in range(numpy.size(qubit_position_j[0, :])):
         if edge_matrix_indices[int(not(qubit_position_j[0, edge_index]))][
                 qubit_position_j[1, edge_index]] < i:
-                    operator += ((int(qubit_position_j[1, edge_index]), 'Z'),)
+            operator += ((int(qubit_position_j[1, edge_index]), 'Z'),)
     a_ij += QubitOperator(operator, 1)
     if j < i:
         a_ij = -1*a_ij
diff --git a/src/openfermion/transforms/_conversion.py b/src/openfermion/transforms/_conversion.py
index 1d51171..662b970 100644
--- a/src/openfermion/transforms/_conversion.py
+++ b/src/openfermion/transforms/_conversion.py
@@ -306,7 +306,7 @@ def get_diagonal_coulomb_hamiltonian(fermion_operator,
                                      n_qubits=None,
                                      ignore_incompatible_terms=False):
     r"""Convert a FermionOperator to a DiagonalCoulombHamiltonian.
-    
+
     Args:
         fermion_operator(FermionOperator): The operator to convert.
         n_qubits(int): Optionally specify the total number of qubits in the
diff --git a/src/openfermion/transforms/_verstraete_cirac.py b/src/openfermion/transforms/_verstraete_cirac.py
index 101a8dc..3b66460 100644
--- a/src/openfermion/transforms/_verstraete_cirac.py
+++ b/src/openfermion/transforms/_verstraete_cirac.py
@@ -17,7 +17,7 @@ import networkx
 
 from openfermion.ops import FermionOperator, QubitOperator
 from openfermion.transforms import jordan_wigner
-from openfermion.utils import majorana_operator 
+from openfermion.utils import majorana_operator
 
 
 def verstraete_cirac_2d_square(operator, x_dimension, y_dimension,
diff --git a/src/openfermion/utils/_lattice.py b/src/openfermion/utils/_lattice.py
index 4227fcc..d512005 100644
--- a/src/openfermion/utils/_lattice.py
+++ b/src/openfermion/utils/_lattice.py
@@ -23,7 +23,7 @@ class Spin(IntEnum):
 
 class SpinPairs(Enum):
     """The spin orbitals corresponding to a pair of spatial orbitals."""
-    
+
     SAME = 1
     ALL = 0
     DIFF = -1
@@ -83,7 +83,7 @@ class HubbardLattice:
     @abc.abstractmethod
     def site_pairs_iter(self, edge_type, ordered=True):
         """Iterable over pairs of sites corresponding to the given edge type."""
-    
+
 
     # properties
 
@@ -113,7 +113,7 @@ class HubbardLattice:
 
 
     def to_spin_orbital_index(self, site_index, dof_index, spin_index):
-        """The index of the spin orbital.""" 
+        """The index of the spin orbital."""
         return (self.site_index_offset(site_index) +
                 self.dof_index_offset(dof_index) +
                 spin_index)
@@ -186,7 +186,7 @@ class HubbardLattice:
 
 class HubbardSquareLattice(HubbardLattice):
     r"""A square lattice for a Hubbard model.
-    
+
     Valid edge types are:
         * 'onsite'
         * 'horizontal_neighbor'
@@ -229,7 +229,7 @@ class HubbardSquareLattice(HubbardLattice):
     @property
     def n_sites(self):
         return self.x_dimension * self.y_dimension
-    
+
 
     @property
     def edge_types(self):
@@ -281,25 +281,23 @@ class HubbardSquareLattice(HubbardLattice):
     def n_horizontal_neighbor_pairs(self, ordered=True):
         """Number of horizontally neighboring (unordered) pairs of sites."""
         n_horizontal_edges_per_y = (
-                self.x_dimension - 
-                (self.x_dimension <= 2 or not self.periodic))
-        return (self.y_dimension * n_horizontal_edges_per_y * 
+            self.x_dimension - (self.x_dimension <= 2 or not self.periodic))
+        return (self.y_dimension * n_horizontal_edges_per_y *
                 (2 if ordered else 1))
 
 
     def n_vertical_neighbor_pairs(self, ordered=True):
         """Number of vertically neighboring (unordered) pairs of sites."""
-        n_vertical_edges_per_x = (
-                self.y_dimension - 
-                (self.y_dimension <= 2 or not self.periodic))
-        return (self.x_dimension * n_vertical_edges_per_x * 
+        n_vertical_edges_per_x = (self.y_dimension -
+                                  (self.y_dimension <= 2 or not self.periodic))
+        return (self.x_dimension * n_vertical_edges_per_x *
                 (2 if ordered else 1))
 
     def n_neighbor_pairs(self, ordered=True):
         """Number of neighboring (unordered) pairs of sites."""
-        return (self.n_horizontal_neighbor_pairs(ordered) + 
+        return (self.n_horizontal_neighbor_pairs(ordered) +
                 self.n_vertical_neighbor_pairs(ordered))
-    
+
 
     def neighbors_iter(self, ordered=True):
         return itertools.chain(
@@ -307,9 +305,8 @@ class HubbardSquareLattice(HubbardLattice):
                 self.vertical_neighbors_iter(ordered))
 
     def diagonal_neighbors_iter(self, ordered=True):
-        n_sites_per_y = (
-                self.x_dimension - 
-                (self.x_dimension <= 2 or not self.periodic))
+        n_sites_per_y = (self.x_dimension -
+                         (self.x_dimension <= 2 or not self.periodic))
         n_sites_per_x = (
                 self.y_dimension -
                 (self.y_dimension <= 2 or not self.periodic))
@@ -325,8 +322,7 @@ class HubbardSquareLattice(HubbardLattice):
 
     def horizontal_neighbors_iter(self, ordered=True):
         n_horizontal_edges_per_y = (
-                self.x_dimension - 
-                (self.x_dimension <= 2 or not self.periodic))
+            self.x_dimension - (self.x_dimension <= 2 or not self.periodic))
         for x in range(n_horizontal_edges_per_y):
             for y in range(self.y_dimension):
                 i = self.to_site_index((x, y))
@@ -337,9 +333,8 @@ class HubbardSquareLattice(HubbardLattice):
 
 
     def vertical_neighbors_iter(self, ordered=True):
-        n_vertical_edges_per_x = (
-                self.y_dimension - 
-                (self.y_dimension <= 2 or not self.periodic))
+        n_vertical_edges_per_x = (self.y_dimension -
+                                  (self.y_dimension <= 2 or not self.periodic))
         for y in range(n_vertical_edges_per_x):
             for x in range(self.x_dimension):
                 i = self.to_site_index((x, y))
diff --git a/src/openfermion/utils/_operator_utils.py b/src/openfermion/utils/_operator_utils.py
index f36bdf5..c62f8e4 100644
--- a/src/openfermion/utils/_operator_utils.py
+++ b/src/openfermion/utils/_operator_utils.py
@@ -971,35 +971,38 @@ def group_into_tensor_product_basis_sets(operator, seed=None):
 
     return sub_operators
 
-  
+
 def _commutes(operator1,operator2):
     return operator1*operator2==operator2*operator1
 
 def _non_fully_commuting_terms(hamiltonian):
     terms = list([QubitOperator(key) for key in hamiltonian.terms.keys()])
-    T=[] # will contain the subset of terms that do not commute universally in terms                                                                                                 
+    T = [
+    ]  # will contain the subset of terms that do not commute universally in terms
     for i in range(len(terms)):
         if any(not _commutes(terms[i],terms[j]) for j in range(len(terms))):
             T.append(terms[i])
     return T
 
 def is_contextual(hamiltonian):
-    """                                                                                                                                                                              
-    Determine whether a hamiltonian (instance of QubitOperator) is contextual,                                                                                                       
-    in the sense of https://arxiv.org/abs/1904.02260.                                                                                                                                
-                                                                                                                                                                                     
-    Args:                                                                                                                                                                            
-        hamiltonian (QubitOperator): the hamiltonian whose                                                                                                                           
-            contextuality is to be evaluated.                                                                                                                                        
-                                                                                                                                                                                     
-    Returns:                                                                                                                                                                         
-        boolean indicating whether hamiltonian is contextual or not.                                                                                                                 
+    """
+    Determine whether a hamiltonian (instance of QubitOperator) is contextual,
+    in the sense of https://arxiv.org/abs/1904.02260.
+
+    Args:
+        hamiltonian (QubitOperator): the hamiltonian whose
+            contextuality is to be evaluated.
+
+    Returns:
+        boolean indicating whether hamiltonian is contextual or not.
     """
     T = _non_fully_commuting_terms(hamiltonian)
-    # Search in T for triples in which exactly one pair anticommutes; if any exist, hamiltonian is contextual.                                                                       
-    for i in range(len(T)): # WLOG, i indexes the operator that (putatively) commutes with both others.                                                                              
+    # Search in T for triples in which exactly one pair anticommutes; if any exist, hamiltonian is contextual.
+    for i in range(
+            len(T)
+    ):  # WLOG, i indexes the operator that (putatively) commutes with both others.
         for j in range(len(T)):
-            for k in range(j+1,len(T)): # Ordering of j, k does not matter.                                                                                                          
+            for k in range(j + 1, len(T)):  # Ordering of j, k does not matter.
                 if i!=j and i!=k and _commutes(T[i],T[j]) and _commutes(T[i],T[k]) and not _commutes(T[j],T[k]):
                     return True
-    return False  
+    return False
diff --git a/src/openfermion/utils/_qubit_tapering_from_stabilizer.py b/src/openfermion/utils/_qubit_tapering_from_stabilizer.py
index d48a6a9..74fba1f 100644
--- a/src/openfermion/utils/_qubit_tapering_from_stabilizer.py
+++ b/src/openfermion/utils/_qubit_tapering_from_stabilizer.py
@@ -1,496 +1,491 @@
-#   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.
-    """
-    if isinstance(stabilizers, (list, tuple, numpy.ndarray)):
-        n_qbits_stabs = 0
-        for ent in stabilizers:
-            if count_qubits(ent) > n_qbits_stabs:
-                n_qbits_stabs = count_qubits(ent)
-    else:
-        n_qbits_stabs = count_qubits(stabilizers)
-
-    n_qbits = max(count_qubits(operator), n_qbits_stabs)
-
-    (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
+#   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.
+    """
+    if isinstance(stabilizers, (list, tuple, numpy.ndarray)):
+        n_qbits_stabs = 0
+        for ent in stabilizers:
+            if count_qubits(ent) > n_qbits_stabs:
+                n_qbits_stabs = count_qubits(ent)
+    else:
+        n_qbits_stabs = count_qubits(stabilizers)
+
+    n_qbits = max(count_qubits(operator), n_qbits_stabs)
+
+    (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
diff --git a/src/openfermion/utils/_sparse_tools.py b/src/openfermion/utils/_sparse_tools.py
index 788a588..35fdd97 100644
--- a/src/openfermion/utils/_sparse_tools.py
+++ b/src/openfermion/utils/_sparse_tools.py
@@ -499,7 +499,7 @@ def jw_sz_restrict_state(state, sz_value,
 
 def jw_get_ground_state_at_particle_number(sparse_operator, particle_number):
     """Compute ground energy and state at a specified particle number.
-    
+
     Assumes the Jordan-Wigner transform. The input operator should be Hermitian
     and particle-number-conserving.
 


Test patch:
diff --git a/src/openfermion/hamiltonians/_general_hubbard_test.py b/src/openfermion/hamiltonians/_general_hubbard_test.py
index cebcfed..eab026a 100644
--- a/src/openfermion/hamiltonians/_general_hubbard_test.py
+++ b/src/openfermion/hamiltonians/_general_hubbard_test.py
@@ -26,26 +26,27 @@ from openfermion.hamiltonians._general_hubbard import InteractionParameter
 def fermi_hubbard_from_general(x_dimension, y_dimension, tunneling, coulomb,
                   chemical_potential=0.,
                   periodic=True, spinless=False, magnetic_field=0):
-    lattice = HubbardSquareLattice(x_dimension, y_dimension, 
-                             periodic=periodic, spinless=spinless)
+    lattice = HubbardSquareLattice(x_dimension,
+                                   y_dimension,
+                                   periodic=periodic,
+                                   spinless=spinless)
     interaction_edge_type = 'neighbor' if spinless else 'onsite'
-    model = FermiHubbardModel(lattice, 
-            tunneling_parameters=(('neighbor', (0, 0), tunneling),),
-            interaction_parameters=((interaction_edge_type, (0, 0), coulomb),),
-            potential_parameters=((0, chemical_potential),),
-            magnetic_field=magnetic_field)
+    model = FermiHubbardModel(lattice,
+                              tunneling_parameters=(('neighbor', (0, 0),
+                                                     tunneling),),
+                              interaction_parameters=((interaction_edge_type,
+                                                       (0, 0), coulomb),),
+                              potential_parameters=((0, chemical_potential),),
+                              magnetic_field=magnetic_field)
     return model.hamiltonian()
 
 
 @pytest.mark.parametrize(
-        'x_dimension,y_dimension,tunneling,coulomb,' + 
-        'chemical_potential,spinless,periodic,magnetic_field',
-        itertools.product(
-            range(1, 4), range(1, 4),
-            (random.uniform(0, 2.),), (random.uniform(0, 2.),),
-            (random.uniform(0, 2.),), (True, False), (True, False),
-            (random.uniform(-1, 1),))
-        )
+    'x_dimension,y_dimension,tunneling,coulomb,' +
+    'chemical_potential,spinless,periodic,magnetic_field',
+    itertools.product(range(1, 4), range(1, 4), (random.uniform(0, 2.),),
+                      (random.uniform(0, 2.),), (random.uniform(0, 2.),),
+                      (True, False), (True, False), (random.uniform(-1, 1),)))
 def test_fermi_hubbard_square_special_general_equivalence(
         x_dimension, y_dimension, tunneling, coulomb,
         chemical_potential, spinless, periodic, magnetic_field):
@@ -61,8 +62,9 @@ def test_fermi_hubbard_square_special_general_equivalence(
 
 def random_parameters(lattice, probability=0.5, distinguish_edges=False):
     parameters = {}
-    edge_types = (('onsite', 'horizontal_neighbor', 'vertical_neighbor') if 
-                  distinguish_edges else ('onsite', 'neighbor'))
+    edge_types = (('onsite', 'horizontal_neighbor',
+                   'vertical_neighbor') if distinguish_edges else
+                  ('onsite', 'neighbor'))
 
     parameters['tunneling_parameters'] = [
             (edge_type, dofs, random.uniform(-1, 1))
@@ -72,12 +74,12 @@ def random_parameters(lattice, probability=0.5, distinguish_edges=False):
 
     possible_spin_pairs = (SpinPairs.ALL,) if lattice.spinless else (SpinPairs.SAME, SpinPairs.DIFF)
     parameters['interaction_parameters'] = [
-            (edge_type, dofs, random.uniform(-1, 1), spin_pairs)
-            for edge_type in edge_types
-            for spin_pairs in possible_spin_pairs
-            for dofs in lattice.dof_pairs_iter(edge_type == 'onsite' and 
-                                               spin_pairs in (SpinPairs.ALL, SpinPairs.SAME))
-            if random.random() <= probability]
+        (edge_type, dofs, random.uniform(-1, 1), spin_pairs)
+        for edge_type in edge_types for spin_pairs in possible_spin_pairs
+        for dofs in lattice.dof_pairs_iter(edge_type == 'onsite' and spin_pairs
+                                           in (SpinPairs.ALL, SpinPairs.SAME))
+        if random.random() <= probability
+    ]
 
     parameters['potential_parameters'] = [
             (dof, random.uniform(-1, 1))
@@ -86,7 +88,7 @@ def random_parameters(lattice, probability=0.5, distinguish_edges=False):
 
     if random.random() <= probability:
         parameters['magnetic_field'] = random.uniform(-1, 1)
-    
+
     return parameters
 
 
@@ -137,12 +139,15 @@ def test_fermi_hubbard_square_lattice_random_parameters(
         if len(spin_orbitals) == 2:
             (i, a, s), (ii, aa, ss) = (
                     lattice.from_spin_orbital_index(i) for i in spin_orbitals)
-            edge_type = ({(0, 0): 'onsite', (0, 1): 'vertical_neighbor', 
-                          (1, 0): 'horizontal_neighbor'}
-                         if distinguish_edges else 
-                         {(0, 0): 'onsite', (0, 1): 'neighbor',
-                          (1, 0): 'neighbor'}
-                        )[lattice.delta_mag(i, ii, True)]
+            edge_type = ({
+                (0, 0): 'onsite',
+                (0, 1): 'vertical_neighbor',
+                (1, 0): 'horizontal_neighbor'
+            } if distinguish_edges else {
+                (0, 0): 'onsite',
+                (0, 1): 'neighbor',
+                (1, 0): 'neighbor'
+            })[lattice.delta_mag(i, ii, True)]
             dofs = tuple(sorted((a, aa)))
             if len(term) == 2:
                 parameter = (edge_type, dofs, -coefficient)
@@ -150,8 +155,8 @@ def test_fermi_hubbard_square_lattice_random_parameters(
                 terms_per_parameter['tunneling', parameter] += 1
             else:
                 assert len(term) == 4
-                spin_pairs = (SpinPairs.ALL if lattice.spinless else 
-                           (SpinPairs.SAME if s == ss else SpinPairs.DIFF))
+                spin_pairs = (SpinPairs.ALL if lattice.spinless else
+                              (SpinPairs.SAME if s == ss else SpinPairs.DIFF))
                 parameter = (edge_type, dofs, coefficient, spin_pairs)
                 assert parameter in parameters['interaction_parameters']
                 terms_per_parameter['interaction', parameter] += 1
@@ -163,19 +168,19 @@ def test_fermi_hubbard_square_lattice_random_parameters(
             potential_coefficient = -coefficient
             if not lattice.spinless:
                 spin = (-1) ** spin_index
-                potential_coefficient -= (
-                        ((-1) ** spin_index) * 
-                         parameters.get('magnetic_field', 0))
+                potential_coefficient -= (((-1)**spin_index) *
+                                          parameters.get('magnetic_field', 0))
             if not potential_coefficient:
                 continue
             parameter = (dof, potential_coefficient)
             assert parameter in parameters['potential_parameters']
             terms_per_parameter['potential', parameter] += 1
     edge_type_to_n_site_pairs = {
-            'onsite': lattice.n_sites, 
-            'neighbor': lattice.n_neighbor_pairs(False),
-            'vertical_neighbor': lattice.n_vertical_neighbor_pairs(False),
-            'horizontal_neighbor': lattice.n_horizontal_neighbor_pairs(False)}
+        'onsite': lattice.n_sites,
+        'neighbor': lattice.n_neighbor_pairs(False),
+        'vertical_neighbor': lattice.n_vertical_neighbor_pairs(False),
+        'horizontal_neighbor': lattice.n_horizontal_neighbor_pairs(False)
+    }
     for (term_type, parameter), n_terms in terms_per_parameter.items():
         if term_type == 'potential':
             assert n_terms == lattice.n_sites * lattice.n_spin_values
@@ -191,8 +196,8 @@ def test_fermi_hubbard_square_lattice_random_parameters(
                 expected_n_terms *= len(set(parameter.dofs))
             if not lattice.spinless:
                 assert parameter.spin_pairs != SpinPairs.ALL
-                if (parameter.edge_type == 'onsite' and 
-                      parameter.spin_pairs == SpinPairs.DIFF):
+                if (parameter.edge_type == 'onsite' and
+                        parameter.spin_pairs == SpinPairs.DIFF):
                     expected_n_terms *= len(set(parameter.dofs))
                 else:
                     expected_n_terms *= 2
diff --git a/src/openfermion/ops/_ising_operator_test.py b/src/openfermion/ops/_ising_operator_test.py
index d14a856..20e50ad 100644
--- a/src/openfermion/ops/_ising_operator_test.py
+++ b/src/openfermion/ops/_ising_operator_test.py
@@ -30,7 +30,6 @@ class GeneralTest(unittest.TestCase):
     def test_ising_operator(self):
         equals_tester = EqualsTester(self)
 
-        group = [IsingOperator('Z0 Z3'), 
-                 IsingOperator([(0, 'Z'), (3, 'Z')])]
+        group = [IsingOperator('Z0 Z3'), IsingOperator([(0, 'Z'), (3, 'Z')])]
 
         equals_tester.add_equality_group(*group)
diff --git a/src/openfermion/utils/_lattice_test.py b/src/openfermion/utils/_lattice_test.py
index c6956b3..a1b43ea 100644
--- a/src/openfermion/utils/_lattice_test.py
+++ b/src/openfermion/utils/_lattice_test.py
@@ -22,15 +22,12 @@ def test_spin():
     lattice = HubbardSquareLattice(3, 3)
     assert tuple(lattice.spin_indices) == (Spin.UP, Spin.DOWN)
 
-@pytest.mark.parametrize(
-        "x_dimension,y_dimension,n_dofs,spinless,periodic",
-        itertools.product(
-            random.sample(range(3, 10), 3), 
-            random.sample(range(3, 10), 3), 
-            range(1, 4),
-            (False, True),
-            (False, True)
-        ))
+
+@pytest.mark.parametrize("x_dimension,y_dimension,n_dofs,spinless,periodic",
+                         itertools.product(random.sample(range(3, 10), 3),
+                                           random.sample(range(3, 10), 3),
+                                           range(1, 4), (False, True),
+                                           (False, True)))
 def test_hubbard_square_lattice(
     x_dimension, y_dimension, n_dofs, spinless, periodic):
     lattice = HubbardSquareLattice(x_dimension, y_dimension,
@@ -41,8 +38,8 @@ def test_hubbard_square_lattice(
     site_indices = tuple(lattice.to_site_index(site) for site in sites)
     assert (sites == tuple(
         lattice.from_site_index(site_index) for site_index in site_indices))
-    assert (site_indices == tuple(lattice.site_indices) == 
-            tuple(range(x_dimension * y_dimension)))
+    assert (site_indices == tuple(lattice.site_indices) == tuple(
+        range(x_dimension * y_dimension)))
 
     spin_orbital_composite_indices = tuple(itertools.product(
         range(x_dimension),range(y_dimension),
@@ -60,9 +57,8 @@ def test_hubbard_square_lattice(
             (x_dimension * (y_dimension - (not periodic))) +
             ((x_dimension - (not periodic)) * y_dimension))
     neighbor_pairs = tuple(lattice.site_pairs_iter('neighbor'))
-    assert (2 * len(tuple(lattice.site_pairs_iter('neighbor', False))) 
-            == len(neighbor_pairs)
-            == n_neighbor_pairs)
+    assert (2 * len(tuple(lattice.site_pairs_iter('neighbor', False))) ==
+            len(neighbor_pairs) == n_neighbor_pairs)
     for i, j in neighbor_pairs:
         assert sum(lattice.delta_mag(i, j, True)) == 1
 
@@ -129,8 +125,8 @@ def test_hubbard_square_lattice_1xd(d):
                 0)
 
 
-@pytest.mark.parametrize('x,y', 
-        (random.sample(range(3, 10), 2) for _ in range(3)))
+@pytest.mark.parametrize('x,y',
+                         (random.sample(range(3, 10), 2) for _ in range(3)))
 def test_hubbard_square_lattice_neighbors(x, y):
     for periodic in (True, False):
         lattice = HubbardSquareLattice(x, y, periodic=periodic)
@@ -155,10 +151,10 @@ def test_hubbard_square_lattice_2xd(d):
                 2 * len(tuple(lattice.neighbors_iter(False))) ==
                 n_neighbor_pairs)
         n_diagonal_neighbor_pairs = 4 * (d - (not periodic))
-        assert (len(tuple(lattice.diagonal_neighbors_iter())) == 
+        assert (len(tuple(lattice.diagonal_neighbors_iter())) ==
                 2 * len(tuple(lattice.diagonal_neighbors_iter(False))) ==
                 n_diagonal_neighbor_pairs)
-               
+
 
 def test_hubbard_square_lattice_2x2():
     for periodic in (True, False):
@@ -191,8 +187,9 @@ def test_spin_pairs_iter():
     with pytest.raises(ValueError):
         spinful_lattice.spin_pairs_iter(10)
 
-    assert (tuple(spinful_lattice.spin_pairs_iter(SpinPairs.ALL, True)) == 
-            ((0, 0), (0, 1), (1, 0), (1, 1)))
+    assert (tuple(spinful_lattice.spin_pairs_iter(SpinPairs.ALL,
+                                                  True)) == ((0, 0), (0, 1),
+                                                             (1, 0), (1, 1)))
     assert (tuple(spinful_lattice.spin_pairs_iter(SpinPairs.ALL, False)) ==
             ((0, 0), (0, 1), (1, 1)))
     assert (tuple(spinful_lattice.spin_pairs_iter(SpinPairs.SAME, True)) ==
@@ -209,7 +206,6 @@ def test_spin_pairs_iter():
             tuple(spinless_lattice.spin_pairs_iter(SpinPairs.SAME, True)) ==
             tuple(spinless_lattice.spin_pairs_iter(SpinPairs.SAME, False)) ==
             ((0, 0),))
-    assert (tuple(spinless_lattice.spin_pairs_iter(SpinPairs.DIFF, True)) ==
-            tuple(spinless_lattice.spin_pairs_iter(SpinPairs.DIFF, False)) ==
-            tuple())
-
+    assert (tuple(spinless_lattice.spin_pairs_iter(
+        SpinPairs.DIFF, True)) == tuple(
+            spinless_lattice.spin_pairs_iter(SpinPairs.DIFF, False)) == tuple())
diff --git a/src/openfermion/utils/_qubit_tapering_from_stabilizer_test.py b/src/openfermion/utils/_qubit_tapering_from_stabilizer_test.py
index 00600e4..4c81ed1 100644
--- a/src/openfermion/utils/_qubit_tapering_from_stabilizer_test.py
+++ b/src/openfermion/utils/_qubit_tapering_from_stabilizer_test.py
@@ -1,350 +1,333 @@
-#   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, count_qubits
-
-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])
-
-    def test_tappering_stabilizer_more_qubits(self):
-        """Test for stabilizer with more qubits than operator."""
-        hamiltonian = QubitOperator('Y0 Y1', 1.0)
-        stab = QubitOperator('X0 X1 X2', -1.0)
-
-        num_qubits = max(count_qubits(hamiltonian),
-                         count_qubits(stab))
-        tap_ham = taper_off_qubits(hamiltonian, stab)
-        num_qubits_tap = count_qubits(tap_ham)
-
-        self.assertFalse(num_qubits == num_qubits_tap)
+#   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, count_qubits
+
+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])
+
+    def test_tappering_stabilizer_more_qubits(self):
+        """Test for stabilizer with more qubits than operator."""
+        hamiltonian = QubitOperator('Y0 Y1', 1.0)
+        stab = QubitOperator('X0 X1 X2', -1.0)
+
+        num_qubits = max(count_qubits(hamiltonian), count_qubits(stab))
+        tap_ham = taper_off_qubits(hamiltonian, stab)
+        num_qubits_tap = count_qubits(tap_ham)
+
+        self.assertFalse(num_qubits == num_qubits_tap)