1""" Hamiltonian formulation of problems """
2from itertools import product
3import numpy as np
4
5from qiskit import QuantumCircuit
6from qiskit.quantum_info import SparsePauliOp
7from qiskit_optimization.converters import QuadraticProgramToQubo
8from qiskit_optimization.translators import from_ising
9
10from qlauncher.base import formatter
11from qlauncher.base.adapter_structure import adapter
12import qlauncher.problems.problem_initialization as problems
13import qlauncher.hampy as hampy
14from qlauncher.hampy import Equation, Variable
15from qlauncher.problems.problem_formulations.hamiltonians.tsp import problem_to_hamiltonian as tsp_to_hamiltonian
16
17
[docs]
18@adapter('hamiltonian', 'qubo', onehot='quadratic')
19def hamiltonian_to_qubo(hamiltonian):
20 qp = from_ising(hamiltonian)
21 conv = QuadraticProgramToQubo()
22 qubo = conv.convert(qp).objective
23 return qubo.quadratic.to_array(), 0
24
25
[docs]
26@adapter("qubo", "hamiltonian")
27def qubo_to_hamiltonian(qubo: np.ndarray) -> SparsePauliOp:
28 q_matrix, offset = qubo
29 num_vars = q_matrix.shape[0]
30 pauli = 0
31 for i, col in enumerate(q_matrix):
32 for j, entry in enumerate(col):
33 if entry == 0:
34 continue
35 if i == j:
36 pauli += SparsePauliOp.from_sparse_list([('I', [0], .5), ('Z', [i], -.5)], num_vars)*entry
37 else:
38 pauli += SparsePauliOp.from_sparse_list([('I', [0], .25), ('Z', [i], -.25),
39 ('Z', [j], -.25), ('ZZ', [i, j], .25)], num_vars)*entry
40 pauli += SparsePauliOp.from_sparse_list([('I', [], offset)], num_vars)
41 return pauli
42
43
[docs]
44def ring_ham(ring: set, n):
45 total = None
46 ring = list(ring)
47 for index in range(len(ring) - 1):
48 sparse_list = []
49 sparse_list.append((("XX", [ring[index], ring[index + 1]], 1)))
50 sparse_list.append((("YY", [ring[index], ring[index + 1]], 1)))
51 sp = SparsePauliOp.from_sparse_list(sparse_list, n)
52 if total is None:
53 total = sp
54 else:
55 total += sp
56 sparse_list = []
57 sparse_list.append((("XX", [ring[-1], ring[0]], 1)))
58 sparse_list.append((("YY", [ring[-1], ring[0]], 1)))
59 sp = SparsePauliOp.from_sparse_list(sparse_list, n)
60 total += sp
61 return SparsePauliOp(total)
62
63
64@formatter(problems.EC, 'hamiltonian')
65class ECQiskit:
66 def __call__(self, problem: problems.EC, onehot='exact') -> SparsePauliOp:
67 """ generating hamiltonian"""
68 elements = set().union(*problem.instance)
69 onehots = []
70 for ele in elements:
71 ohs = set()
72 for i, subset in enumerate(problem.instance):
73 if ele in subset:
74 ohs.add(i)
75 onehots.append(ohs)
76 hamiltonian = None
77 for ohs in onehots:
78 if onehot == 'exact':
79 part = (~hampy.one_in_n(list(ohs), len(problem.instance))).hamiltonian
80 elif onehot == 'quadratic':
81 part = hampy.one_in_n(list(ohs), len(problem.instance), quadratic=True).hamiltonian
82
83 if hamiltonian is None:
84 hamiltonian = part
85 else:
86 hamiltonian += part
87 return hamiltonian.simplify()
88
89 def get_mixer_hamiltonian(self, problem: problems.EC, amount_of_rings=None):
90 """ generates mixer hamiltonian """
91 def get_main_set():
92 main_set = []
93 for element_set in problem.instance:
94 for elem in element_set:
95 if elem not in main_set:
96 main_set.append(elem)
97 return main_set
98
99 def get_constraints():
100 constraints, main_set = [], get_main_set()
101 for element in main_set:
102 element_set = set()
103 for index, _ in enumerate(problem.instance):
104 if element in problem.instance[index]:
105 element_set.add(index)
106 if len(element_set) > 0 and element_set not in constraints:
107 constraints.append(element_set)
108
109 return constraints
110
111 # creating mixer hamiltonians for all qubits that aren't in rings (in other words applying X gate to them)
112 def x_gate_ham(x_gate: list):
113 total = None
114 for elem in x_gate:
115 sparse_list = []
116 sparse_list.append((("X", [elem], 1)))
117 sp = SparsePauliOp.from_sparse_list(
118 sparse_list, len(problem.instance))
119 if total is None:
120 total = sp
121 else:
122 total += sp
123 return SparsePauliOp(total)
124
125 # looking for all rings in a data and creating a list with them
126 ring, x_gate, constraints = [], [], get_constraints()
127
128 ring.append(max(constraints, key=len))
129
130 ring_qubits = set.union(*ring)
131
132 for set_ in constraints:
133 if len(set_.intersection(ring_qubits)) == 0:
134 ring.append(set_)
135 ring_qubits.update(set_)
136
137 if amount_of_rings is not None:
138 max_amount_of_rings, user_rings = len(ring), []
139 if amount_of_rings > max_amount_of_rings:
140 raise ValueError(
141 f"Too many rings. Maximum amount is {max_amount_of_rings}")
142 elif amount_of_rings == 0:
143 ring_qubits = []
144 else:
145 current_qubits = ring[0]
146 for index in range(amount_of_rings):
147 user_rings.append(ring[index])
148 current_qubits = current_qubits.union(ring[index])
149 ring_qubits = current_qubits
150 x_gate.extend(id for id, _ in enumerate(
151 problem.instance) if id not in ring_qubits)
152
153 # connecting all parts of mixer hamiltonian together
154 mix_ham = None
155 for set_ in ring:
156 if mix_ham is None:
157 mix_ham = ring_ham(set_, len(problem.instance))
158 else:
159 mix_ham += ring_ham(set_, len(problem.instance))
160
161 if mix_ham is None:
162 mix_ham = x_gate_ham(x_gate)
163 else:
164 mix_ham += x_gate_ham(x_gate)
165
166 return mix_ham
167
168
169@formatter(problems.JSSP, 'hamiltonian')
170def get_qiskit_hamiltonian(problem: problems.JSSP) -> SparsePauliOp:
171 if problem.optimization_problem:
172 return problem.h_o
173 else:
174 return problem.h_d
175
176
177@formatter(problems.MaxCut, 'hamiltonian')
178def get_qiskit_hamiltonian(problem: problems.MaxCut):
179 ham = None
180 n = problem.instance.number_of_nodes()
181 for edge in problem.instance.edges():
182 if ham is None:
183 ham = ~hampy.one_in_n(edge, n)
184 else:
185 ham += ~hampy.one_in_n(edge, n)
186 return ham.hamiltonian.simplify()
187
188
189@formatter(problems.QATM, 'hamiltonian')
190class QATMQiskit:
191 def __call__(self, problem: problems.QATM, onehot='exact') -> SparsePauliOp:
192 cm = problem.instance['cm']
193 aircrafts = problem.instance['aircrafts']
194
195 onehot_hamiltonian = None
196 for plane, manouvers in aircrafts.groupby(by='aircraft'):
197 if onehot == 'exact':
198 h = (~hampy.one_in_n(manouvers.index.values.tolist(), len(cm))).hamiltonian
199 elif onehot == 'quadratic':
200 h = hampy.one_in_n(manouvers.index.values.tolist(), len(cm), quadratic=True).hamiltonian
201 elif onehot == 'xor':
202 total = None
203 eq = Equation(len(cm))
204 for part in manouvers.index.values.tolist():
205 if total is None:
206 total = eq[part].to_equation()
207 continue
208 total ^= eq[part]
209 h = (~total).hamiltonian
210 if onehot_hamiltonian is not None:
211 onehot_hamiltonian += h
212 else:
213 onehot_hamiltonian = h
214
215 triu = np.triu(cm, k=1)
216 conflict_hamiltonian = None
217 for p1, p2 in zip(*np.where(triu == 1)):
218 eq = Equation(len(cm))
219 partial_hamiltonian = (eq[int(p1)] & eq[int(p2)]).hamiltonian
220 if conflict_hamiltonian is not None:
221 conflict_hamiltonian += partial_hamiltonian
222 else:
223 conflict_hamiltonian = partial_hamiltonian
224
225 hamiltonian = onehot_hamiltonian + conflict_hamiltonian
226
227 if problem.optimization_problem:
228 goal_hamiltonian = None
229 for i, (maneuver, ac) in problem.instance['aircrafts'].iterrows():
230 if maneuver != ac:
231 eq = Equation(len(aircrafts))
232 h = Variable(i, eq).to_equation()
233 if goal_hamiltonian is None:
234 goal_hamiltonian = h
235 else:
236 goal_hamiltonian += h
237 goal_hamiltonian /= sum(sum(cm))
238 hamiltonian += goal_hamiltonian
239
240 return hamiltonian.simplify()
241
242 def get_mixer_hamiltonian(self, problem: problems.QATM) -> SparsePauliOp:
243 cm = problem.instance['cm']
244 aircrafts = problem.instance['aircrafts']
245
246 mixer_hamiltonian = None
247 for plane, manouvers in aircrafts.groupby(by='aircraft'):
248 h = ring_ham(manouvers.index.values.tolist(), len(cm))
249 if mixer_hamiltonian is None:
250 mixer_hamiltonian = h
251 else:
252 mixer_hamiltonian += h
253 return mixer_hamiltonian
254
255 def get_QAOAAnsatz_initial_state(self, problem: problems.QATM) -> QuantumCircuit:
256 aircrafts = problem.instance['aircrafts']
257 qc = QuantumCircuit(len(aircrafts))
258 for plane, manouvers in aircrafts.groupby(by='aircraft'):
259 qc.x(manouvers.index.values.tolist()[0])
260 return qc
261
262
263@formatter(problems.Raw, 'hamiltonian')
264def get_qiskit_hamiltonian(problem: problems.Raw) -> SparsePauliOp:
265 return problem.instance
266
267
268@formatter(problems.TSP, 'hamiltonian')
269def get_qiskit_hamiltonian(problem: problems.TSP, onehot='exact', constraints_weight=1, costs_weight=1) -> SparsePauliOp:
270 return tsp_to_hamiltonian(
271 problem,
272 onehot=onehot,
273 constraints_weight=constraints_weight,
274 costs_weight=costs_weight
275 )
276
277
[docs]
278@formatter(problems.GraphColoring, 'hamiltonian')
279def get_qiskit_hamiltonian(problem: problems.GraphColoring, constraints_weight=1, costs_weight=1):
280 color_bit_length = int(np.ceil(np.log2(problem.num_colors)))
281 num_qubits = problem.instance.number_of_nodes() * color_bit_length
282 eq = Equation(num_qubits)
283 # Penalty for assigning the same colors to neighboring vertices
284 for node1, node2 in problem.instance.edges:
285 for ind, comb in enumerate(product(range(2), repeat=color_bit_length)):
286 if ind >= problem.num_colors:
287 break
288 eq2 = None
289 for i in range(color_bit_length):
290 qubit1 = eq[node1 * color_bit_length + i]
291 qubit2 = eq[node2 * color_bit_length + i]
292 if comb[i]:
293 exp = qubit1 & qubit2
294 else:
295 exp = ~qubit1 & ~qubit2
296 if eq2 is None:
297 eq2 = exp
298 else:
299 eq2 &= exp
300 eq += eq2
301 eq *= costs_weight
302 # Penalty for using excessive colors
303 for node in problem.instance.nodes:
304 for ind, comb in enumerate(product(range(2), repeat=color_bit_length)):
305 if ind < problem.num_colors:
306 continue
307 eq2 = None
308 for i in range(color_bit_length):
309 qubit = eq[node * color_bit_length + i]
310 exp = qubit if comb[i] else ~qubit
311 if eq2 is None:
312 eq2 = exp
313 else:
314 eq2 &= exp
315 eq += eq2
316 eq *= constraints_weight
317 return eq.hamiltonian
