Repository: quantumlib/OpenFermion
Pull request #509 (2019-03-19 04:14:37)

adds InteractionOperator utilities
handled by normal_ordered, is_hermitian, and hermitian_conjugate

new zero construction method

new projection method

bryano: @kevinsung I think this is good to go
bryano: @kevinsung Could you merge this?

Patch:
diff --git a/src/openfermion/ops/_interaction_operator.py b/src/openfermion/ops/_interaction_operator.py
index 76a5681..d1be1aa 100644
--- a/src/openfermion/ops/_interaction_operator.py
+++ b/src/openfermion/ops/_interaction_operator.py
@@ -13,6 +13,8 @@
 """Class and functions to store interaction operators."""
 import itertools
 
+import numpy
+
 from openfermion.ops import PolynomialTensor
 
 
@@ -117,6 +119,20 @@ class InteractionOperator(PolynomialTensor):
                 seen |= set(_symmetric_two_body_terms(quad, complex_valued))
                 yield tuple(zip(quad, (1, 1, 0, 0)))
 
+    @classmethod
+    def zero(cls, n_qubits):
+        return cls(0, numpy.zeros((n_qubits,) * 2), numpy.zeros((n_qubits,) * 4))
+
+    def projected(self, indices, exact=False):
+        projected_n_body_tensors = self.projected_n_body_tensors(
+                indices, exact)
+        return type(self)(*(projected_n_body_tensors[key]
+            for key in [(), (1, 0), (1, 1, 0, 0)]))
+
+    def with_function_applied_elementwise(self, func):
+        return type(self)(*(func(tensor) for tensor in
+            [self.constant, self.one_body_tensor, self.two_body_tensor]))
+
 
 def _symmetric_two_body_terms(quad, complex_valued):
     p, q, r, s = quad
diff --git a/src/openfermion/ops/_polynomial_tensor.py b/src/openfermion/ops/_polynomial_tensor.py
index 3eb7273..ea926b2 100644
--- a/src/openfermion/ops/_polynomial_tensor.py
+++ b/src/openfermion/ops/_polynomial_tensor.py
@@ -15,6 +15,8 @@ from __future__ import division
 
 import copy
 import itertools
+import operator
+
 import numpy
 
 from openfermion.config import EQ_TOLERANCE
@@ -221,11 +223,17 @@ class PolynomialTensor(object):
         summand += addend
         return summand
 
-    def __neg__(self):
-        neg_n_body_tensors = dict()
+    def with_function_applied_elementwise(self, func):
+        new_n_body_tensors = dict()
         for key in self.n_body_tensors:
-            neg_n_body_tensors[key] = numpy.negative(self.n_body_tensors[key])
-        return PolynomialTensor(neg_n_body_tensors)
+            new_n_body_tensors[key] = func(self.n_body_tensors[key])
+        return PolynomialTensor(new_n_body_tensors)
+
+    def __neg__(self):
+        return self.with_function_applied_elementwise(operator.neg)
+
+    def __mod__(self, other):
+        return self.with_function_applied_elementwise(lambda x: x % other)
 
     def __isub__(self, subtrahend):
         if not issubclass(type(subtrahend), PolynomialTensor):
@@ -346,3 +354,35 @@ class PolynomialTensor(object):
 
     def __repr__(self):
         return str(self)
+
+    def projected_n_body_tensors(self, selection, exact=False):
+        """Keep only selected elements.
+
+        Args:
+            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. 
+            exact (bool): Whether or not the selection is strict.
+        """
+        comparator = (operator.eq if exact else operator.le)
+        if isinstance(selection, int):
+            pred = lambda index: comparator(len(set(index)), selection)
+            dims = range(self.n_qubits)
+        else:
+            selection = set(selection)
+            pred = lambda index: comparator(set(index), selection)
+            dims = selection
+
+        projected_n_body_tensors = dict()
+        for key, tensor in self.n_body_tensors.items():
+            if not key:
+                projected_n_body_tensors[key] = (
+                        tensor if not (exact and selection) else 0)
+                continue
+            projected_tensor = numpy.zeros_like(tensor)
+            for index in itertools.product(dims, repeat=len(key)):
+                if pred(index):
+                    projected_tensor[index] = tensor[index]
+            projected_n_body_tensors[key] = projected_tensor
+        return projected_n_body_tensors
diff --git a/src/openfermion/utils/_operator_utils.py b/src/openfermion/utils/_operator_utils.py
index 2eb242a..bb8b212 100644
--- a/src/openfermion/utils/_operator_utils.py
+++ b/src/openfermion/utils/_operator_utils.py
@@ -13,17 +13,16 @@
 """This module provides generic tools for classes in ops/"""
 from __future__ import absolute_import
 from builtins import map, zip
-
-import os
 import copy
+import itertools
 import marshal
-from numpy.random import RandomState
+import os
 
 import numpy
+from numpy.random import RandomState
 from scipy.sparse import spmatrix
 
 from openfermion.config import DATA_DIRECTORY, EQ_TOLERANCE
-
 from openfermion.ops import (BosonOperator,
                              DiagonalCoulombHamiltonian,
                              FermionOperator,
@@ -238,6 +237,16 @@ def hermitian_conjugated(operator):
                 sorted(conjugate_term, key=lambda factor: factor[0]))
             conjugate_operator.terms[conjugate_term] = coefficient.conjugate()
 
+    # Handle InteractionOperator
+    elif isinstance(operator, InteractionOperator):
+        conjugate_constant = operator.constant.conjugate()
+        conjugate_one_body_tensor = hermitian_conjugated(
+            operator.one_body_tensor)
+        conjugate_two_body_tensor = hermitian_conjugated(
+            operator.two_body_tensor)
+        conjugate_operator = type(operator)(conjugate_constant,
+            conjugate_one_body_tensor, conjugate_two_body_tensor)
+
     # Handle sparse matrix
     elif isinstance(operator, spmatrix):
         conjugate_operator = operator.getH()
@@ -256,8 +265,9 @@ def hermitian_conjugated(operator):
 
 def is_hermitian(operator):
     """Test if operator is Hermitian."""
-    # Handle FermionOperator and BosonOperator
-    if isinstance(operator, (FermionOperator, BosonOperator)):
+    # Handle FermionOperator, BosonOperator, and InteractionOperator
+    if isinstance(operator,
+            (FermionOperator, BosonOperator, InteractionOperator)):
         return (normal_ordered(operator) ==
                 normal_ordered(hermitian_conjugated(operator)))
 
@@ -800,7 +810,7 @@ def normal_ordered_quad_term(term, coefficient, hbar=1.):
 
 def normal_ordered(operator, hbar=1.):
     r"""Compute and return the normal ordered form of a FermionOperator,
-    BosonOperator, QuadOperator.
+    BosonOperator, QuadOperator, or InteractionOperator.
 
     Due to the canonical commutation/anticommutation relations satisfied
     by these operators, there are multiple forms that the same operator
@@ -817,7 +827,7 @@ def normal_ordered(operator, hbar=1.):
 
     Args:
         operator: an instance of the FermionOperator, BosonOperator,
-            or QuadOperator classes.
+            QuadOperator, or InteractionOperator classes.
         hbar (float): the value of hbar used in the definition of the
             commutator [q_i, p_j] = i hbar delta_ij. By default hbar=1.
             This argument only applies when normal ordering QuadOperators.
@@ -839,9 +849,36 @@ def normal_ordered(operator, hbar=1.):
         order_fn = normal_ordered_quad_term
         kwargs['hbar'] = hbar
 
+    elif isinstance(operator, InteractionOperator):
+        constant = operator.constant
+        n_modes = operator.n_qubits
+        one_body_tensor = operator.one_body_tensor.copy()
+        two_body_tensor = numpy.zeros_like(operator.two_body_tensor)
+        quadratic_index_pairs = (
+            (pq, pq) for pq in itertools.combinations(range(n_modes)[::-1], 2))
+        cubic_index_pairs = (index_pair
+            for p, q, r in itertools.combinations(range(n_modes)[::-1], 3)
+            for index_pair in [
+                ((p, q), (p, r)), ((p, r), (p, q)),
+                ((p, q), (q, r)), ((q, r), (p, q)),
+                ((p, r), (q, r)), ((q, r), (p, r))])
+        quartic_index_pairs = (index_pair
+            for p, q, r, s in itertools.combinations(range(n_modes)[::-1], 4)
+            for index_pair in [
+                ((p, q), (r, s)), ((r, s), (p, q)),
+                ((p, r), (q, s)), ((q, s), (p, r)),
+                ((p, s), (q, r)), ((q, r), (p, s))])
+        index_pairs = itertools.chain(
+            quadratic_index_pairs, cubic_index_pairs, quartic_index_pairs)
+        for pq, rs in index_pairs:
+            two_body_tensor[pq + rs] = sum(
+                s * ss * operator.two_body_tensor[pq[::s] + rs[::ss]]
+                for s, ss in itertools.product([-1, 1], repeat=2))
+        return InteractionOperator(constant, one_body_tensor, two_body_tensor)
+
     else:
         raise TypeError('Can only normal order FermionOperator, '
-                        'BosonOperator, or QuadOperator.')
+                        'BosonOperator, QuadOperator, or InteractionOperator.')
 
     for term, coefficient in operator.terms.items():
         ordered_operator += order_fn(term, coefficient, **kwargs)


Test patch:
diff --git a/src/openfermion/ops/_interaction_operator_test.py b/src/openfermion/ops/_interaction_operator_test.py
index 5699757..ff9cf6b 100644
--- a/src/openfermion/ops/_interaction_operator_test.py
+++ b/src/openfermion/ops/_interaction_operator_test.py
@@ -13,6 +13,8 @@
 """Tests for interaction_operator.py."""
 import unittest
 
+import itertools
+
 import numpy
 
 from openfermion.ops import InteractionOperator
@@ -83,3 +85,43 @@ class InteractionOperatorTest(unittest.TestCase):
                                           summed_op.n_body_tensors[(1, 0)]))
         self.assertTrue(numpy.array_equal(summed_op.two_body_tensor,
                                           summed_op.n_body_tensors[(1, 1, 0, 0)]))
+
+    def test_neg(self):
+        interaction_op = InteractionOperator(
+                0, numpy.ones((3, 3)), numpy.ones((3, 3, 3, 3)))
+        neg_interaction_op = -interaction_op
+        assert isinstance(neg_interaction_op, InteractionOperator)
+        assert neg_interaction_op == InteractionOperator(
+                0, -numpy.ones((3, 3)), -numpy.ones((3, 3, 3, 3)))
+
+    def test_zero(self):
+        interaction_op = InteractionOperator(0, numpy.ones((3, 3)), numpy.ones((3, 3, 3, 3)))
+        assert InteractionOperator.zero(interaction_op.n_qubits) == interaction_op * 0
+
+    def test_projected(self):
+        interaction_op = InteractionOperator(1, numpy.ones((3, 3)), numpy.ones((3, 3, 3, 3)))
+        projected_two_body_tensor = numpy.zeros_like(interaction_op.two_body_tensor)
+        for i in range(3):
+            projected_two_body_tensor[(i,) * 4] = 1
+        projected_op = InteractionOperator(0, numpy.eye(3), projected_two_body_tensor)
+        assert projected_op == interaction_op.projected(1, exact=True)
+        projected_op.constant = 1
+        assert projected_op == interaction_op.projected(1, exact=False)
+
+        projected_op.one_body_tensor = numpy.zeros_like(interaction_op.one_body_tensor)
+        for pq in itertools.product(range(2), repeat=2):
+            projected_op.one_body_tensor[pq] = 1
+        projected_op.two_body_tensor = numpy.zeros_like(interaction_op.two_body_tensor)
+        for pqrs in itertools.product(range(2), repeat=4):
+            projected_op.two_body_tensor[pqrs] = 1
+        assert projected_op == interaction_op.projected((0, 1), exact=False)
+
+        projected_op.constant = 0
+        projected_op.one_body_tensor = numpy.zeros_like(interaction_op.one_body_tensor)
+        projected_op.one_body_tensor[0, 1] = 1
+        projected_op.one_body_tensor[1, 0] = 1
+        projected_op.two_body_tensor = numpy.zeros_like(interaction_op.two_body_tensor)
+        for pqrs in itertools.product(range(2), repeat=4):
+            if len(set(pqrs)) > 1:
+                projected_op.two_body_tensor[pqrs] = 1
+        assert projected_op == interaction_op.projected((0, 1), exact=True)
diff --git a/src/openfermion/ops/_polynomial_tensor_test.py b/src/openfermion/ops/_polynomial_tensor_test.py
index 662af92..4e83a5b 100644
--- a/src/openfermion/ops/_polynomial_tensor_test.py
+++ b/src/openfermion/ops/_polynomial_tensor_test.py
@@ -205,6 +205,17 @@ class PolynomialTensorTest(unittest.TestCase):
         new_tensor = self.polynomial_tensor_a + self.polynomial_tensor_b
         self.assertEqual(new_tensor, self.polynomial_tensor_ab)
 
+    def test_mod(self):
+        new_constant = 2.0
+        new_one_body = numpy.zeros_like(self.one_body_a)
+        new_one_body[0, 1] = 2
+        new_two_body = numpy.zeros_like(self.two_body_a)
+        new_two_body[0, 1, 0, 1] = 1
+        new_two_body[1, 1, 0, 0] = 2
+        new_tensor = PolynomialTensor({(): new_constant,
+            (1, 0): new_one_body, (1, 1, 0, 0): new_two_body})
+        assert new_tensor == (self.polynomial_tensor_a % 3)
+
     def test_iadd(self):
         new_tensor = copy.deepcopy(self.polynomial_tensor_a)
         new_tensor += self.polynomial_tensor_b
diff --git a/src/openfermion/utils/_operator_utils_test.py b/src/openfermion/utils/_operator_utils_test.py
index 98562fd..dfcca0b 100644
--- a/src/openfermion/utils/_operator_utils_test.py
+++ b/src/openfermion/utils/_operator_utils_test.py
@@ -13,6 +13,7 @@
 """Tests for operator_utils."""
 from __future__ import absolute_import
 
+import itertools
 import os
 import unittest
 
@@ -30,6 +31,7 @@ from openfermion.utils import (Grid, is_hermitian,
                                random_interaction_operator)
 
 from openfermion.utils._operator_utils import *
+from openfermion.utils._testing_utils import random_interaction_operator
 
 
 class OperatorUtilsTest(unittest.TestCase):
@@ -403,6 +405,14 @@ class HermitianConjugatedTest(unittest.TestCase):
                  BosonOperator('2^ 2', -0.1j))
         self.assertEqual(op_hc, hermitian_conjugated(op))
 
+    def test_hermitian_conjugated_interaction_operator(self):
+        for n_orbitals, _ in itertools.product((1, 2, 5), range(5)):
+            operator = random_interaction_operator(n_orbitals)
+            qubit_operator = jordan_wigner(operator)
+            conjugate_operator = hermitian_conjugated(operator)
+            conjugate_qubit_operator = jordan_wigner(conjugate_operator)
+            assert hermitian_conjugated(qubit_operator) == conjugate_qubit_operator
+
     def test_exceptions(self):
         with self.assertRaises(TypeError):
             _ = is_hermitian('a')
@@ -840,6 +850,33 @@ class TestNormalOrdering(unittest.TestCase):
         self.assertTrue(op_132 == normal_ordered(op_132))
         self.assertTrue(op_132 == normal_ordered(op_321))
 
+    def test_interaction_operator(self):
+        for n_orbitals, real, _ in itertools.product(
+                (1, 2, 5), (True, False), range(5)):
+            operator = random_interaction_operator(n_orbitals, real=real)
+            normal_ordered_operator = normal_ordered(operator)
+            expected_qubit_operator = jordan_wigner(operator)
+            actual_qubit_operator = jordan_wigner(
+                    normal_ordered_operator)
+            assert expected_qubit_operator == actual_qubit_operator
+            two_body_tensor = normal_ordered_operator.two_body_tensor
+            n_orbitals = len(two_body_tensor)
+            ones = numpy.ones((n_orbitals,) * 2)
+            triu = numpy.triu(ones, 1)
+            shape = (n_orbitals ** 2, 1)
+            mask = (triu.reshape(shape) * ones.reshape(shape[::-1]) +
+                    ones.reshape(shape) * triu.reshape(shape[::-1])
+                    ).reshape((n_orbitals,) * 4)
+            assert numpy.allclose(mask * two_body_tensor,
+                    numpy.zeros((n_orbitals,) * 4))
+            for term in normal_ordered_operator:
+                order = len(term) // 2
+                left_term, right_term = term[:order], term[order:]
+                assert all(i[1] == 1 for i in left_term)
+                assert all(i[1] == 0 for i in right_term)
+                assert left_term == tuple(sorted(left_term, reverse=True))
+                assert right_term == tuple(sorted(right_term, reverse=True))
+
     def test_exceptions(self):
         with self.assertRaises(TypeError):
             _ = normal_ordered(1)