Repository: quantumlib/OpenFermion (python)

Issue:
Create class to represent Majorana operators
The algebra of Majorana operators is simple because they either commute or anticommute. Algebraically, they are easier to manipulate than the fermionic creation and annihilation operators; the square of any single-mode Majorana operator always equals one. Any two different single-mode Majorana operators always anticommute.

Pull request description:
Define MajoranaOperator
Fixes #458 

Patch:
diff --git a/src/openfermion/ops/__init__.py b/src/openfermion/ops/__init__.py
index d604a08..202a65b 100644
--- a/src/openfermion/ops/__init__.py
+++ b/src/openfermion/ops/__init__.py
@@ -13,6 +13,7 @@
 from ._binary_polynomial import BinaryPolynomial
 from ._diagonal_coulomb_hamiltonian import DiagonalCoulombHamiltonian
 from ._indexing import down_index, up_index
+from ._majorana_operator import MajoranaOperator
 from ._polynomial_tensor import PolynomialTensor, general_basis_change
 from ._quadratic_hamiltonian import QuadraticHamiltonian
 from ._symbolic_operator import SymbolicOperator
diff --git a/src/openfermion/ops/_majorana_operator.py b/src/openfermion/ops/_majorana_operator.py
new file mode 100644
index 0000000..a4488fc
--- /dev/null
+++ b/src/openfermion/ops/_majorana_operator.py
@@ -0,0 +1,328 @@
+#   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.
+
+"""The MajoranaOperator data structure."""
+
+from __future__ import division
+from future.utils import viewkeys
+
+import numpy
+
+
+class MajoranaOperator:
+    """A linear combination of products of Majorana operators.
+
+    A system of N fermionic modes can be described using 2N Majorana operators
+    :math:`\gamma_1, \ldots, \gamma_{2N}`
+    as an alternative to using N fermionic annihilation operators. The algebra
+    of Majorana operators amounts to the relation
+
+    .. math::
+        \{\gamma_i, \gamma_j\}
+         = \gamma_i \gamma_j + \gamma_j \gamma_i
+         = 2 \delta_{ij}
+
+    Note that this implies :math:`\gamma_i^2 = 1`.
+
+    The MajoranaOperator class stores a linear combination of products
+    of Majorana operators. Each product is represented as a tuple of
+    integers representing the indices of the operators. As an example,
+    `MajoranaOperator((2, 3, 5), -1.5)` initializes an operator with
+    a single term which represents the operator
+    :math:`-1.5 \gamma_2 \gamma_3 \gamma_5`. MajoranaOperators can be
+    added, subtracted, multiplied, and divided by scalars. They can be
+    compared for approximate numerical equality using `==`.
+
+    Attributes:
+        terms: A dictionary from term, represented by a tuple of integers,
+        to the coefficient of the term in the linear combination.
+    """
+
+    def __init__(self, term=None, coefficient=1.0):
+        """Initialize a MajoranaOperator with a single term.
+
+        Args:
+            term (Tuple[int]): The indices of a Majorana operator term
+                to start off with
+            coefficient (complex): The coefficient of the term
+
+        Returns:
+            MajoranaOperator
+        """
+        self.terms = {}
+        if term is not None:
+            term, parity = _sort_majorana_term(term)
+            self.terms[term] = coefficient * (-1)**parity
+
+    @staticmethod
+    def from_dict(terms):
+        """Initialize a MajoranaOperator from a terms dictionary.
+
+        WARNING: The given dictionary is not validated whatsoever. It's up
+        to you to ensure that it is properly formed.
+
+        Args:
+            terms: A dictionary from Majorana term to coefficient
+        """
+        op = MajoranaOperator()
+        op.terms = terms
+        return op
+
+    def commutes_with(self, other):
+        """Test commutation with another MajoranaOperator"""
+        if not isinstance(other, type(self)):
+            raise TypeError(
+                'Can only test commutation with another MajoranaOperator.')
+
+        if len(self.terms) == 1 and len(other.terms) == 1:
+            return _majorana_terms_commute(list(self.terms.keys())[0],
+                                           list(other.terms.keys())[0])
+        return self*other == other*self
+
+    def __iadd__(self, other):
+        if not isinstance(other, type(self)):
+            return NotImplemented
+
+        for term, coefficient in other.terms.items():
+            if term in self.terms:
+                self.terms[term] += coefficient
+            else:
+                self.terms[term] = coefficient
+
+        return self
+
+    def __add__(self, other):
+        if not isinstance(other, type(self)):
+            return NotImplemented
+
+        terms = {}
+        terms.update(self.terms)
+
+        for term, coefficient in other.terms.items():
+            if term in terms:
+                terms[term] += coefficient
+            else:
+                terms[term] = coefficient
+
+        return MajoranaOperator.from_dict(terms)
+
+    def __isub__(self, other):
+        if not isinstance(other, type(self)):
+            return NotImplemented
+
+        for term, coefficient in other.terms.items():
+            if term in self.terms:
+                self.terms[term] -= coefficient
+            else:
+                self.terms[term] = coefficient
+
+        return self
+
+    def __sub__(self, other):
+        if not isinstance(other, type(self)):
+            return NotImplemented
+
+        terms = {}
+        terms.update(self.terms)
+        for term, coefficient in other.terms.items():
+            if term in terms:
+                terms[term] -= coefficient
+            else:
+                terms[term] = -coefficient
+        return MajoranaOperator.from_dict(terms)
+
+    def __mul__(self, other):
+        if not isinstance(other, (type(self), int, float, complex)):
+            return NotImplemented
+
+        if isinstance(other, (int, float, complex)):
+            terms = {term: coefficient*other
+                     for term, coefficient in self.terms.items()}
+            return MajoranaOperator.from_dict(terms)
+
+        terms = {}
+        for left_term, left_coefficient in self.terms.items():
+            for right_term, right_coefficient in other.terms.items():
+                new_term, parity = _merge_majorana_terms(left_term, right_term)
+                coefficient = left_coefficient*right_coefficient*(-1)**parity
+                if new_term in terms:
+                    terms[new_term] += coefficient
+                else:
+                    terms[new_term] = coefficient
+        return MajoranaOperator.from_dict(terms)
+
+    def __imul__(self, other):
+        if not isinstance(other, (type(self), int, float, complex)):
+            return NotImplemented
+
+        if isinstance(other, (int, float, complex)):
+            for term in self.terms:
+                self.terms[term] *= other
+            return self
+
+        return self * other
+
+    def __rmul__(self, other):
+        if not isinstance(other, (int, float, complex)):
+            return NotImplemented
+        return self * other
+
+    def __truediv__(self, other):
+        if not isinstance(other, (int, float, complex)):
+            return NotImplemented
+
+        terms = {term: coefficient/other
+                 for term, coefficient in self.terms.items()}
+        return MajoranaOperator.from_dict(terms)
+
+    def __itruediv__(self, other):
+        if not isinstance(other, (int, float, complex)):
+            return NotImplemented
+
+        for term in self.terms:
+            self.terms[term] /= other
+        return self
+
+    def __div__(self, divisor):
+        """ For compatibility with Python 2. """
+        return self.__truediv__(divisor)
+
+    def __pow__(self, other):
+        if not isinstance(other, int):
+            return NotImplemented
+
+        if other < 0:
+            raise TypeError('Cannot raise MajoranaOperator to negative power.')
+
+        result = MajoranaOperator(())
+        for _ in range(other):
+            result *= self
+        return result
+
+    def __neg__(self):
+        return -1 * self
+
+    def __eq__(self, other):
+        """Approximate numerical equality."""
+        if not isinstance(other, type(self)):
+            return NotImplemented
+
+        for term in viewkeys(self.terms) | viewkeys(other.terms):
+            if term in self.terms and term in other.terms:
+                if not numpy.isclose(self.terms[term], other.terms[term]):
+                    return False
+            elif term in self.terms:
+                if not numpy.isclose(self.terms[term], 0.0):
+                    return False
+            else:
+                if not numpy.isclose(other.terms[term], 0.0):
+                    return False
+        return True
+
+    def __ne__(self, other):
+        return not self == other
+
+    def __str__(self):
+        if not self.terms:
+            return '0'
+        lines = []
+        for term, coeff in sorted(self.terms.items()):
+            if numpy.isclose(coeff, 0.0):
+                continue
+            lines.append('{} {} +'.format(coeff, term))
+        if not lines:
+            return '0'
+        return '\n'.join(lines)[:-2]
+
+    def __repr__(self):
+        return 'MajoranaOperator.from_dict(terms={!r})'.format(self.terms)
+
+
+def _sort_majorana_term(term):
+    """Sort a Majorana term.
+
+    Args:
+        term (Tuple[int]): The indices of a Majorana operator term
+
+    Returns:
+        Tuple[Tuple[int], int]. The first object returned is a sorted list
+        representing the indices acted upon. The second object is the parity
+        of the term. A parity of 1 indicates that the term should include
+        a minus sign.
+    """
+    if len(term) < 2:
+        return term, 0
+    center = len(term) // 2
+    left_term, left_parity = _sort_majorana_term(term[:center])
+    right_term, right_parity = _sort_majorana_term(term[center:])
+    merged_term, merge_parity = _merge_majorana_terms(left_term, right_term)
+    return merged_term, (left_parity + right_parity + merge_parity) % 2
+
+
+def _merge_majorana_terms(left_term, right_term):
+    """Merge two Majorana terms.
+
+    Args:
+        left_term (Tuple[int]): The left-hand term
+        right_term (Tuple[int]): The right-hand term
+
+    Returns:
+        Tuple[Tuple[int], int]. The first object returned is a sorted list
+        representing the indices acted upon. The second object is the parity
+        of the term. A parity of 1 indicates that the term should include
+        a minus sign.
+    """
+    merged_term = []
+    parity = 0
+    i, j = 0, 0
+    while i < len(left_term) and j < len(right_term):
+        if left_term[i] < right_term[j]:
+            merged_term.append(left_term[i])
+            i += 1
+        elif left_term[i] > right_term[j]:
+            merged_term.append(right_term[j])
+            j += 1
+            parity += len(left_term) - i
+        else:
+            parity += len(left_term) - i - 1
+            i += 1
+            j += 1
+    if i == len(left_term):
+        merged_term.extend(right_term[j:])
+    else:
+        merged_term.extend(left_term[i:])
+    return tuple(merged_term), parity % 2
+
+
+def _majorana_terms_commute(term_a, term_b):
+    """Whether two Majorana terms commute.
+
+    Args:
+        term_a (Tuple[int]): The indices of a Majorana operator term
+        term_b (Tuple[int]): The indices of a Majorana operator term
+
+    Returns:
+        bool. Whether The terms commute.
+    """
+    intersection = 0
+    i, j = 0, 0
+    while i < len(term_a) and j < len(term_b):
+        if term_a[i] < term_b[j]:
+            i += 1
+        elif term_a[i] > term_b[j]:
+            j += 1
+        else:
+            intersection += 1
+            i += 1
+            j += 1
+    parity = (len(term_a)*len(term_b) - intersection) % 2
+    return not parity


Test patch:
diff --git a/src/openfermion/ops/_majorana_operator_test.py b/src/openfermion/ops/_majorana_operator_test.py
new file mode 100644
index 0000000..37d35f8
--- /dev/null
+++ b/src/openfermion/ops/_majorana_operator_test.py
@@ -0,0 +1,200 @@
+#   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.
+
+import pytest
+
+from openfermion import MajoranaOperator
+
+
+def test_majorana_operator_init():
+    op = MajoranaOperator((0, 7, 4, 1, 10))
+    assert op.terms == {(0, 1, 4, 7, 10): -1.0}
+
+    op = MajoranaOperator((0, 1, 0, 1, 0))
+    assert op.terms == {(0,): -1.0}
+
+    op = MajoranaOperator((3, 2, 5, 1, 3, 4))
+    assert op.terms == {(1, 2, 4, 5): 1.0}
+
+    op = MajoranaOperator((5, 10, 4, 3, 6, 9, 6))
+    assert op.terms == {(3, 4, 5, 9, 10): -1.0}
+
+
+def test_majorana_operator_commutes_with():
+    a = MajoranaOperator((0, 1, 5))
+    b = MajoranaOperator((1, 2, 7))
+    c = MajoranaOperator((2, 3, 4))
+    d = MajoranaOperator((0, 3, 6))
+
+    assert a.commutes_with(b)
+    assert not a.commutes_with(c)
+    assert a.commutes_with(d)
+    assert b.commutes_with(c)
+    assert not b.commutes_with(d)
+    assert c.commutes_with(d)
+    assert (a+c).commutes_with(b+d)
+
+    e = MajoranaOperator((0, 1, 1, 1, 4, 5))
+    f = MajoranaOperator((0, 1, 1, 4))
+
+    assert e.commutes_with(f)
+
+    with pytest.raises(TypeError):
+        _ = e.commutes_with(0)
+
+
+
+def test_majorana_operator_add_subtract():
+    a = MajoranaOperator((0, 2, 3), -1.25)
+    b = MajoranaOperator((0, 2, 3), -0.5)
+    c = MajoranaOperator((1, 5, 7), 4.75)
+    d = MajoranaOperator((3, 5, 7), 2.25)
+
+    a += c
+    assert a.terms == {(0, 2, 3): -1.25,
+                       (1, 5, 7): 4.75}
+
+    a -= d
+    assert a.terms == {(0, 2, 3): -1.25,
+                       (1, 5, 7): 4.75,
+                       (3, 5, 7): 2.25}
+
+    a += MajoranaOperator((0, 2, 3), 0.5)
+    assert a.terms == {(0, 2, 3): -.75,
+                       (1, 5, 7): 4.75,
+                       (3, 5, 7): 2.25}
+
+    a -= MajoranaOperator((0, 2, 3), 0.25)
+    assert a.terms == {(0, 2, 3): -1.0,
+                       (1, 5, 7): 4.75,
+                       (3, 5, 7): 2.25}
+
+    assert (a + b).terms == {(0, 2, 3): -1.5,
+                             (1, 5, 7): 4.75,
+                             (3, 5, 7): 2.25}
+
+    assert (a - b).terms == {(0, 2, 3): -0.5,
+                             (1, 5, 7): 4.75,
+                             (3, 5, 7): 2.25}
+
+    with pytest.raises(TypeError):
+        _ = a + 0
+
+    with pytest.raises(TypeError):
+        a += 0
+
+    with pytest.raises(TypeError):
+        _ = a - 0
+
+    with pytest.raises(TypeError):
+        a -= 0
+
+
+def test_majorana_operator_multiply():
+    a = MajoranaOperator((0, 1, 5), 1.5) + MajoranaOperator((1, 2, 7), -0.5)
+    b = MajoranaOperator((2, 3, 4), 1.75) - MajoranaOperator((0, 3, 6), 0.25)
+    assert (a * a).terms == {(): -2.5,
+                             (0, 2, 5, 7): -1.5}
+    assert (a * b).terms == {(1, 3, 4, 7): 0.875,
+                             (0, 1, 2, 3, 6, 7): -0.125,
+                             (1, 3, 5, 6): 0.375,
+                             (0, 1, 2, 3, 4, 5): -2.625}
+    assert (2 * a).terms == (a * 2).terms == {(0, 1, 5): 3.0,
+                                              (1, 2, 7): -1.0}
+
+    a *= 2
+    a *= MajoranaOperator(())
+    assert a.terms == {(0, 1, 5): 3.0,
+                       (1, 2, 7): -1.0}
+
+    with pytest.raises(TypeError):
+        _ = a * 'a'
+
+    with pytest.raises(TypeError):
+        _ = 'a' * a
+
+    with pytest.raises(TypeError):
+        a *= 'a'
+
+
+def test_majorana_operator_pow():
+    a = MajoranaOperator((0, 1, 5), 1.5) + MajoranaOperator((1, 2, 7), -0.5)
+    assert (a**2).terms == {(): -2.5,
+                            (0, 2, 5, 7): -1.5}
+
+    with pytest.raises(TypeError):
+        _ = a**-1
+
+    with pytest.raises(TypeError):
+        _ = a**'a'
+
+
+def test_majorana_operator_divide():
+    a = MajoranaOperator((0, 1, 5), 1.5) + MajoranaOperator((1, 2, 7), -0.5)
+    assert (a / 2).terms == {(0, 1, 5): 0.75,
+                             (1, 2, 7): -0.25}
+
+    a /= 2
+    assert a.terms == {(0, 1, 5): 0.75,
+                       (1, 2, 7): -0.25}
+
+    with pytest.raises(TypeError):
+        _ = a / 'a'
+
+    with pytest.raises(TypeError):
+        a /= 'a'
+
+
+def test_majorana_operator_neg():
+    a = MajoranaOperator((0, 1, 5), 1.5) + MajoranaOperator((1, 2, 7), -0.5)
+    assert (-a).terms == {(0, 1, 5): -1.5,
+                          (1, 2, 7): 0.5}
+
+
+def test_majorana_operator_eq():
+    a = MajoranaOperator((0, 1, 5), 1.5) + MajoranaOperator((1, 2, 7), -0.5)
+    b = (MajoranaOperator((0, 1, 5), 1.5) +
+         MajoranaOperator((1, 2, 7), -0.5) +
+         MajoranaOperator((3, 4, 5), 0.0))
+    c = (MajoranaOperator((0, 1, 5), 1.5) +
+         MajoranaOperator((1, 2, 7), -0.5) +
+         MajoranaOperator((3, 4, 5), 0.1))
+    d = MajoranaOperator((0, 1, 5), 1.75) + MajoranaOperator((1, 2, 7), -0.75)
+    e = MajoranaOperator((0, 1, 5), 1.5) - MajoranaOperator((0, 3, 6), 0.25)
+
+    assert a == b
+    assert a != c
+    assert a != d
+    assert a != e
+
+    assert a != 0
+
+
+def test_majorana_operator_str():
+    zero = MajoranaOperator()
+    one = MajoranaOperator(())
+    still_zero = MajoranaOperator((0, 1, 5), 0.0)
+    a = MajoranaOperator((0, 1, 5), 1.5)
+    b = MajoranaOperator((1, 2, 7), -0.5)
+
+    assert str(zero) == '0'
+    assert str(one) == '1.0 ()'
+    assert str(still_zero) == '0'
+    assert str(a) == '1.5 (0, 1, 5)'
+    assert str(b) == '-0.5 (1, 2, 7)'
+    assert str(a+b) == """1.5 (0, 1, 5) +
+-0.5 (1, 2, 7)"""
+
+
+def test_majorana_operator_repr():
+    a = MajoranaOperator((0, 1, 5), 1.5)
+    assert repr(a) == 'MajoranaOperator.from_dict(terms={(0, 1, 5): 1.5})'