1#   Licensed under the Apache License, Version 2.0 (the "License");
2#   you may not use this file except in compliance with the License.
3#   You may obtain a copy of the License at
4#
5#       http://www.apache.org/licenses/LICENSE-2.0
6#
7#   Unless required by applicable law or agreed to in writing, software
8#   distributed under the License is distributed on an "AS IS" BASIS,
9#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
10#   See the License for the specific language governing permissions and
11#   limitations under the License.
12""" The transform function that does Fermion-qubit mappings
13    based on a BinaryCode (arXiv:1712.07067) """
14
15from typing import Any, List, Tuple, Union
16
17import numpy
18
19from openfermion.ops.operators import (BinaryCode, FermionOperator,
20                                       QubitOperator, BinaryPolynomial)
21
22
23def extractor(binary_op: BinaryPolynomial) -> QubitOperator:
24    """ Applies the extraction superoperator to a binary expression
25     to obtain the corresponding qubit operators.
26
27    Args:
28        binary_op (BinaryPolynomial): the binary term
29
30    Returns (QubitOperator): the qubit operator corresponding to the
31        binary terms
32    """
33    return_fn = 1
34    for term in binary_op.terms:
35        multiplier = 1
36        if len(term) == 1:
37            term = term[0]
38            if isinstance(term, (numpy.int32, numpy.int64, int)):
39                multiplier = QubitOperator('Z' + str(term))
40            else:
41                multiplier = -1
42        elif len(term) > 1:
43            multiplier = dissolve(term)
44
45        return_fn *= multiplier
46    return return_fn
47
48
49def dissolve(term: Tuple[Any, Any]) -> QubitOperator:
50    """Decomposition helper. Takes a product of binary variables
51    and outputs the Pauli-string sum that corresponds to the
52    decomposed multi-qubit operator.
53
54    Args:
55        term (tuple): product of binary variables, i.e.: 'w0 w2 w3'
56
57    Returns (QubitOperator): superposition of Pauli-strings
58
59    Raises:
60        ValueError: if the variable in term is not integer
61    """
62    prod = 2.0
63    for var in term:
64        if not isinstance(var, (numpy.int32, numpy.int64, int)):
65            raise ValueError('dissolve only works on integers')
66        prod *= (QubitOperator((), 0.5) - QubitOperator('Z' + str(var), 0.5))
67    return QubitOperator((), 1.0) - prod
68
69
70def make_parity_list(code: BinaryCode) -> List[BinaryPolynomial]:
71    """Create the parity list from the decoder of the input code.
72    The output parity list has a similar structure as code.decoder.
73
74    Args:
75        code (BinaryCode): the code to extract the parity list from.
76
77    Returns (list): list of BinaryPolynomial, the parity list
78
79    Raises:
80        TypeError: argument is not BinaryCode
81    """
82    if not isinstance(code, BinaryCode):
83        raise TypeError('argument is not BinaryCode')
84    parity_binaries = [BinaryPolynomial()]
85    for index in numpy.arange(code.n_modes - 1):
86        parity_binaries += [parity_binaries[-1] + code.decoder[index]]
87    return parity_binaries
88
89
90def binary_code_transform(hamiltonian: FermionOperator,
91                          code: BinaryCode) -> QubitOperator:
92    """ Transforms a Hamiltonian written in fermionic basis into a Hamiltonian
93    written in qubit basis, via a binary code.
94
95    The role of the binary code is to relate the occupation vectors (v0 v1 v2
96    ... vN-1) that span the fermionic basis, to the qubit basis, spanned by
97    binary vectors (w0, w1, w2, ..., wn-1).
98
99    The binary code has to provide an analytic relation between the binary
100    vectors (v0, v1, ..., vN-1) and (w0, w1, ..., wn-1), and possibly has the
101    property N>n, when the Fermion basis is smaller than the fermionic Fock
102    space. The binary_code_transform function can transform Fermion operators
103    to qubit operators for custom- and qubit-saving mappings.
104
105    Note:
106        Logic multi-qubit operators are decomposed into Pauli-strings (e.g.
107        CPhase(1,2) = 0.5 * (1 + Z1 + Z2 - Z1 Z2 ) ), which might increase
108        the number of Hamiltonian terms drastically.
109
110    Args:
111        hamiltonian (FermionOperator): the fermionic Hamiltonian
112        code (BinaryCode): the binary code to transform the Hamiltonian
113
114    Returns (QubitOperator): the transformed Hamiltonian
115
116    Raises:
117        TypeError: if the hamiltonian is not a FermionOperator or code is not
118        a BinaryCode
119    """
120    if not isinstance(hamiltonian, FermionOperator):
121        raise TypeError('hamiltonian provided must be a FermionOperator'
122                        'received {}'.format(type(hamiltonian)))
123
124    if not isinstance(code, BinaryCode):
125        raise TypeError('code provided must be a BinaryCode'
126                        'received {}'.format(type(code)))
127
128    new_hamiltonian = QubitOperator()
129    parity_list = make_parity_list(code)
130
131    # for each term in hamiltonian
132    for term, term_coefficient in hamiltonian.terms.items():
133        """ the updated parity and occupation account for sign changes due
134        changed occupations mid-way in the term """
135        updated_parity = 0  # parity sign exponent
136        parity_term = BinaryPolynomial()
137        changed_occupation_vector = [0] * code.n_modes
138        transformed_term = QubitOperator(())
139
140        # keep track of indices appeared before
141        fermionic_indices = numpy.array([])
142
143        # for each multiplier
144        for op_idx, op_tuple in enumerate(reversed(term)):
145            # get count exponent, parity exponent addition
146            fermionic_indices = numpy.append(fermionic_indices, op_tuple[0])
147            count = numpy.count_nonzero(
148                fermionic_indices[:op_idx] == op_tuple[0])
149            updated_parity += numpy.count_nonzero(
150                fermionic_indices[:op_idx] < op_tuple[0])
151
152            # update term
153            extracted = extractor(code.decoder[op_tuple[0]])
154            extracted *= (((-1)**count) * ((-1)**(op_tuple[1])) * 0.5)
155            transformed_term *= QubitOperator((), 0.5) - extracted
156
157            # update parity term and occupation vector
158            changed_occupation_vector[op_tuple[0]] += 1
159            parity_term += parity_list[op_tuple[0]]
160
161        # the parity sign and parity term
162        transformed_term *= QubitOperator((), (-1)**updated_parity)
163        transformed_term *= extractor(parity_term)
164
165        # the update operator
166        changed_qubit_vector = numpy.mod(
167            code.encoder.dot(changed_occupation_vector), 2)
168
169        update_operator = QubitOperator(())
170        for index, q_vec in enumerate(changed_qubit_vector):
171            if q_vec:
172                update_operator *= QubitOperator('X' + str(index))
173
174        # append new term to new hamiltonian
175        new_hamiltonian += term_coefficient * update_operator * \
176                           transformed_term
177
178    new_hamiltonian.compress()
179    return new_hamiltonian
180