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