Skip to content
Permalink
Browse files

adds InteractionOperator utilities (#509)

* adds InteractionOperator utilities

handled by normal_ordered, is_hermitian, and hermitian_conjugate

new zero construction method

new projection method

* handle dicts

* reversed ordering for consistency

* added comment for hermitian conjugate of dict

* preserve type of InteractionOperator

* better preservation of type

* refactored

* bug fix

* coverage

* removed handling of dicts

* removed '_and_dict' from test name
  • Loading branch information...
bryano authored and kevinsung committed Apr 29, 2019
1 parent 1b3154e commit a6264653c29f760abe9e99b15c00e25c4936b61f
@@ -13,6 +13,8 @@
"""Class and functions to store interaction operators."""
import itertools

import numpy

from openfermion.ops import PolynomialTensor


@@ -117,6 +119,20 @@ def unique_iter(self, complex_valued=False):
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
@@ -13,6 +13,8 @@
"""Tests for interaction_operator.py."""
import unittest

import itertools

import numpy

from openfermion.ops import InteractionOperator
@@ -83,3 +85,43 @@ def test_addition(self):
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)
@@ -15,6 +15,8 @@

import copy
import itertools
import operator

import numpy

from openfermion.config import EQ_TOLERANCE
@@ -221,11 +223,17 @@ def __add__(self, addend):
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 @@ def rotate_basis(self, rotation_matrix):

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
@@ -205,6 +205,17 @@ def test_add(self):
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
@@ -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)
@@ -13,6 +13,7 @@
"""Tests for operator_utils."""
from __future__ import absolute_import

import itertools
import os
import unittest

@@ -30,6 +31,7 @@
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 @@ def test_hermitian_conjugated_semihermitian(self):
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 @@ def test_quad_triple(self):
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)

0 comments on commit a626465

Please sign in to comment.
You can’t perform that action at this time.