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.
12from typing import Union
13import itertools
14import numpy
15import sympy
16
17from openfermion.ops.operators import (QuadOperator, BosonOperator,
18                                       FermionOperator, MajoranaOperator)
19from openfermion.ops.representations import (PolynomialTensor,
20                                             DiagonalCoulombHamiltonian)
21
22from openfermion.utils.operator_utils import count_qubits
23
24
25def get_quad_operator(operator, hbar=1.):
26    """Convert to QuadOperator.
27
28    Args:
29        operator: BosonOperator.
30        hbar (float): the value of hbar used in the definition
31            of the commutator [q_i, p_j] = i hbar delta_ij.
32            By default hbar=1.
33
34    Returns:
35        quad_operator: An instance of the QuadOperator class.
36    """
37    quad_operator = QuadOperator()
38
39    if isinstance(operator, BosonOperator):
40        for term, coefficient in operator.terms.items():
41            tmp = QuadOperator('', coefficient)
42            for i, d in term:
43                tmp *= (1./numpy.sqrt(2.*hbar)) \
44                    * (QuadOperator(((i, 'q')))
45                        + QuadOperator(((i, 'p')), 1j*(-1)**d))
46            quad_operator += tmp
47
48    else:
49        raise TypeError("Only BosonOperator is currently "
50                        "supported for get_quad_operator.")
51
52    return quad_operator
53
54
55def check_no_sympy(operator):
56    """Checks whether a SymbolicOperator contains any
57    sympy expressions, which will prevent it being converted
58    to a PolynomialTensor or DiagonalCoulombHamiltonian
59
60    Args:
61        operator(SymbolicOperator): the operator to be tested
62    """
63    for key in operator.terms:
64        if isinstance(operator.terms[key], sympy.Expr):
65            raise TypeError('This conversion is currently not supported ' +
66                            'for operators with sympy expressions ' +
67                            'as coefficients')
68
69
70def get_boson_operator(operator, hbar=1.):
71    """Convert to BosonOperator.
72
73    Args:
74        operator: QuadOperator.
75        hbar (float): the value of hbar used in the definition
76            of the commutator [q_i, p_j] = i hbar delta_ij.
77            By default hbar=1.
78
79    Returns:
80        boson_operator: An instance of the BosonOperator class.
81    """
82    boson_operator = BosonOperator()
83
84    if isinstance(operator, QuadOperator):
85        for term, coefficient in operator.terms.items():
86            tmp = BosonOperator('', coefficient)
87            for i, d in term:
88                if d == 'q':
89                    coeff = numpy.sqrt(hbar / 2)
90                    sign = 1
91                elif d == 'p':
92                    coeff = -1j * numpy.sqrt(hbar / 2)
93                    sign = -1
94
95                tmp *= coeff * (BosonOperator(((i, 0))) + BosonOperator(
96                    ((i, 1)), sign))
97            boson_operator += tmp
98
99    else:
100        raise TypeError("Only QuadOperator is currently "
101                        "supported for get_boson_operator.")
102
103    return boson_operator
104
105
106def get_fermion_operator(operator):
107    """Convert to FermionOperator.
108
109    Returns:
110        fermion_operator: An instance of the FermionOperator class.
111    """
112    if isinstance(operator, PolynomialTensor):
113        return _polynomial_tensor_to_fermion_operator(operator)
114    elif isinstance(operator, DiagonalCoulombHamiltonian):
115        return _diagonal_coulomb_hamiltonian_to_fermion_operator(operator)
116    elif isinstance(operator, MajoranaOperator):
117        return _majorana_operator_to_fermion_operator(operator)
118    else:
119        raise TypeError('{} cannot be converted to FermionOperator'.format(
120            type(operator)))
121
122
123def _polynomial_tensor_to_fermion_operator(operator):
124    fermion_operator = FermionOperator()
125    for term in operator:
126        fermion_operator += FermionOperator(term, operator[term])
127    return fermion_operator
128
129
130def _diagonal_coulomb_hamiltonian_to_fermion_operator(operator):
131    fermion_operator = FermionOperator()
132    n_qubits = count_qubits(operator)
133    fermion_operator += FermionOperator((), operator.constant)
134    for p, q in itertools.product(range(n_qubits), repeat=2):
135        fermion_operator += FermionOperator(((p, 1), (q, 0)),
136                                            operator.one_body[p, q])
137        fermion_operator += FermionOperator(((p, 1), (p, 0), (q, 1), (q, 0)),
138                                            operator.two_body[p, q])
139    return fermion_operator
140
141
142def _majorana_operator_to_fermion_operator(majorana_operator):
143    fermion_operator = FermionOperator()
144    for term, coeff in majorana_operator.terms.items():
145        converted_term = _majorana_term_to_fermion_operator(term)
146        converted_term *= coeff
147        fermion_operator += converted_term
148    return fermion_operator
149
150
151def _majorana_term_to_fermion_operator(term):
152    converted_term = FermionOperator(())
153    for index in term:
154        j, b = divmod(index, 2)
155        if b:
156            converted_op = FermionOperator((j, 0), -1j)
157            converted_op += FermionOperator((j, 1), 1j)
158        else:
159            converted_op = FermionOperator((j, 0))
160            converted_op += FermionOperator((j, 1))
161        converted_term *= converted_op
162    return converted_term
163
164
165def get_majorana_operator(
166        operator: Union[PolynomialTensor, DiagonalCoulombHamiltonian,
167                        FermionOperator]) -> MajoranaOperator:
168    """
169    Convert to MajoranaOperator.
170
171    Uses the convention of even + odd indexing of Majorana modes derived from
172    a fermionic mode:
173        fermion annhil.  c_k  -> ( gamma_{2k} + 1.j * gamma_{2k+1} ) / 2
174        fermion creation c^_k -> ( gamma_{2k} - 1.j * gamma_{2k+1} ) / 2
175
176    Args:
177        operator (PolynomialTensor,
178            DiagonalCoulombHamiltonian or
179            FermionOperator): Operator to write as Majorana Operator.
180
181    Returns:
182        majorana_operator: An instance of the MajoranaOperator class.
183
184    Raises:
185        TypeError: If operator is not of PolynomialTensor,
186            DiagonalCoulombHamiltonian or FermionOperator.
187    """
188    if isinstance(operator, FermionOperator):
189        return _fermion_operator_to_majorana_operator(operator)
190    elif isinstance(operator, (PolynomialTensor, DiagonalCoulombHamiltonian)):
191        return _fermion_operator_to_majorana_operator(
192            get_fermion_operator(operator))
193    raise TypeError('{} cannot be converted to MajoranaOperator'.format(
194        type(operator)))
195
196
197def _fermion_operator_to_majorana_operator(fermion_operator: FermionOperator
198                                          ) -> MajoranaOperator:
199    """
200    Convert FermionOperator to MajoranaOperator.
201
202    Auxiliar function of get_majorana_operator.
203
204    Args:
205        fermion_operator (FermionOperator): To convert to MajoranaOperator.
206
207    Returns:
208        majorana_operator object.
209
210    Raises:
211        TypeError: if input is not a FermionOperator.
212    """
213    if not isinstance(fermion_operator, FermionOperator):
214        raise TypeError('Input a FermionOperator.')
215
216    majorana_operator = MajoranaOperator()
217    for term, coeff in fermion_operator.terms.items():
218        converted_term = _fermion_term_to_majorana_operator(term)
219        converted_term *= coeff
220        majorana_operator += converted_term
221
222    return majorana_operator
223
224
225def _fermion_term_to_majorana_operator(term: tuple) -> MajoranaOperator:
226    """
227    Convert single terms of FermionOperator to Majorana.
228    (Auxiliary function of get_majorana_operator.)
229
230    Convention: even + odd indexing of Majorana modes derived from a
231    fermionic mode:
232        fermion annhil.  c_k  -> ( gamma_{2k} + 1.j * gamma_{2k+1} ) / 2
233        fermion creation c^_k -> ( gamma_{2k} - 1.j * gamma_{2k+1} ) / 2
234
235    Args:
236        term (tuple): single FermionOperator term.
237
238    Returns:
239        converted_term: single MajoranaOperator term.
240
241    Raises:
242        TypeError: if term is a tuple.
243    """
244    if not isinstance(term, tuple):
245        raise TypeError('Term does not have the correct Type.')
246
247    converted_term = MajoranaOperator(())
248    for index, action in term:
249        converted_op = MajoranaOperator((2 * index,), 0.5)
250
251        if action:
252            converted_op += MajoranaOperator((2 * index + 1,), -0.5j)
253
254        else:
255            converted_op += MajoranaOperator((2 * index + 1,), 0.5j)
256
257        converted_term *= converted_op
258
259    return converted_term
260