1from collections.abc import Callable
2from typing import TYPE_CHECKING, Any, overload
3
4import numpy as np
5from dimod import BinaryQuadraticModel
6from pyqubo import Binary # type: ignore
7from pyqubo import Model as pyquboModel
8from qiskit import QuantumCircuit, QuantumRegister
9from qiskit.circuit import Gate
10from qiskit.circuit.library import QFTGate, SwapGate, UnitaryGate, XGate
11from qiskit.quantum_info import SparsePauliOp
12from qiskit_nature.second_q.drivers import PySCFDriver
13from qiskit_nature.second_q.formats.molecule_info import MoleculeInfo
14from qiskit_nature.second_q.mappers import ParityMapper, QubitMapper
15from qiskit_optimization.converters import QuadraticProgramToQubo
16from qiskit_optimization.translators import from_ising
17
18if TYPE_CHECKING:
19 from qiskit_nature.second_q.problems import ElectronicStructureProblem
20
21from qlauncher.hampy import Equation
22
23
[docs]
24class Model:
25 _all_problems: dict[str, type['Model']] = {}
26
27 def __init__(self, instance: Any) -> None:
28 self.instance = instance
29
30 def __init_subclass__(cls) -> None:
31 if Model not in cls.__bases__:
32 return
33 Model._all_problems[cls.__name__] = cls
34 cls._mapping: dict[type[Model], Callable[[], Model]] = {}
35 for method_name in cls.__dict__:
36 if method_name.startswith('to_'):
37 method = cls.__dict__[method_name]
38 cls._mapping[method.__annotations__['return']] = method
39
40
[docs]
41class QUBO(Model):
42 def __init__(self, matrix: np.ndarray, offset: float = 0) -> None:
43 self.matrix = (matrix + matrix.T) / 2
44 self.offset = offset
45
[docs]
46 def to_hamiltonian(self) -> 'Hamiltonian':
47 num_vars = self.matrix.shape[0]
48 pauli = 0
49 for i, col in enumerate(self.matrix):
50 for j, entry in enumerate(col):
51 if entry == 0:
52 continue
53 if i == j:
54 pauli += SparsePauliOp.from_sparse_list([('I', [0], 0.5), ('Z', [i], -0.5)], num_vars) * entry
55 else:
56 pauli += (
57 SparsePauliOp.from_sparse_list(
58 [('I', [0], 0.25), ('Z', [i], -0.25), ('Z', [j], -0.25), ('ZZ', [i, j], 0.25)], num_vars
59 )
60 * entry
61 )
62 pauli += SparsePauliOp.from_sparse_list([('I', [0], self.offset)], num_vars)
63 return Hamiltonian(Equation(pauli))
64
[docs]
65 def to_fn(self) -> 'FN':
66 def function(bin_vec: np.ndarray) -> float:
67 return np.dot(bin_vec, np.dot(self.matrix, bin_vec)) + self.offset
68
69 return FN(function)
70
[docs]
71 def to_bqm(self: 'QUBO') -> 'BQM':
72 matrix = self.matrix
73 symmetric = (self.matrix.transpose() == self.matrix).all()
74 if not symmetric:
75 for i in range(len(self.matrix)):
76 for j in range(len(self.matrix)):
77 if i > j:
78 matrix[i][j] = 0
79 elif j > j:
80 matrix[i][j] *= 2
81
82 values_and_qubits = {(x, y): c for y, r in enumerate(matrix) for x, c in enumerate(r) if c != 0}
83 number_of_qubits = len(matrix)
84 qubits = [Binary(f'x{i}') for i in range(number_of_qubits)]
85 H = 0
86 for (x, y), value in values_and_qubits.items():
87 H += value * qubits[x] * qubits[y]
88 model = H.compile()
89 return BQM(model)
90
91
[docs]
92class FN(Model):
93 def __init__(self, function: Callable[[np.ndarray], float]) -> None:
94 self.function = function
95
96 def __call__(self, vector: np.ndarray) -> float:
97 return self.function(vector)
98
99
[docs]
100class Hamiltonian(Model):
101 @overload
102 def __init__(
103 self, hamiltonian: SparsePauliOp, mixer_hamiltonian: SparsePauliOp | None = None, initial_state: QuantumCircuit | None = None
104 ) -> None: ...
105 @overload
106 def __init__(
107 self,
108 hamiltonian: Equation,
109 mixer_hamiltonian: Equation | None = None,
110 initial_state: QuantumCircuit | None = None,
111 ) -> None: ...
112 def __init__(
113 self,
114 hamiltonian: Equation | SparsePauliOp,
115 mixer_hamiltonian: Equation | SparsePauliOp | None = None,
116 initial_state: QuantumCircuit | None = None,
117 ) -> None:
118 if isinstance(hamiltonian, SparsePauliOp):
119 hamiltonian = Equation(hamiltonian)
120 if isinstance(mixer_hamiltonian, SparsePauliOp):
121 mixer_hamiltonian = Equation(mixer_hamiltonian)
122 self._hampy_equation = hamiltonian
123 self._hampy_mixer_equation: Equation | None = mixer_hamiltonian
124 self._initial_state: QuantumCircuit | None = initial_state
125
126 @property
127 def mixer_hamiltonian(self) -> SparsePauliOp | None:
128 if self._hampy_mixer_equation is not None:
129 return self._hampy_mixer_equation.hamiltonian
130 return None
131
132 @property
133 def hamiltonian(self) -> SparsePauliOp:
134 return self._hampy_equation.hamiltonian
135
136 @property
137 def is_quadratic(self) -> bool:
138 return self._hampy_equation.is_quadratic()
139
140 @mixer_hamiltonian.setter
141 def mixer_hamiltonian(self, mixer_hamiltonian: SparsePauliOp) -> None:
142 self._hampy_mixer_equation = Equation(mixer_hamiltonian)
143
144 @property
145 def initial_state(self) -> QuantumCircuit | None:
146 return self._initial_state
147
148 @initial_state.setter
149 def initial_state(self, initial_state: QuantumCircuit) -> None:
150 self._initial_state = initial_state
151
[docs]
152 def to_qubo(self) -> QUBO:
153 qp = from_ising(self.hamiltonian)
154 conv = QuadraticProgramToQubo()
155 qubo = conv.convert(qp).objective
156 return QUBO(qubo.quadratic.to_array(), float(qubo.constant))
157
158
[docs]
159class BQM(Model):
160 def __init__(self, model: pyquboModel | BinaryQuadraticModel) -> None: # noqa: ANN401
161 if isinstance(model, pyquboModel):
162 self.model = model
163 self._bqm = None
164 else:
165 self.model = None
166 self._bqm = model
167
168 @property
169 def bqm(self) -> BinaryQuadraticModel:
170 if self._bqm is not None:
171 return self._bqm
172 self._bqm = self.model.to_bqm()
173 return self._bqm
174
[docs]
175 def to_qubo(self) -> QUBO:
176 """Returns Qubo function"""
177 if self.model is None:
178 raise NotImplementedError('To transfer into QUBO it is required to provide to BQM pyqubo model')
179 model = self.model
180 variables = sorted(model.variables)
181 num_qubits = len(variables)
182 Q_matrix = np.zeros((num_qubits, num_qubits))
183 inv_map = {v: i for i, v in enumerate(variables)}
184 qubo_dict, offset = model.to_qubo()
185 for key, value in qubo_dict.items():
186 var1, var2 = key
187 Q_matrix[inv_map[var1], inv_map[var2]] = value
188 return QUBO(Q_matrix, offset)
189
[docs]
190 def to_hamiltonian(self) -> Hamiltonian:
191 """Returns Hamiltonian function"""
192
193 bqm = self.model.to_bqm()
194 variables, new_offset = bqm.variables, bqm.offset
195 variables = list(variables)
196 variables.sort()
197
198 sparse_list = []
199
200 for i, coeff in bqm.linear.items():
201 sparse_list.append(('Z', [variables.index(i)], -coeff / 2))
202 new_offset += coeff / 2
203
204 for (i, j), coeff in bqm.quadratic.items():
205 sparse_list.append(('ZZ', [variables.index(i), variables.index(j)], coeff / 4))
206 sparse_list.append(('Z', [variables.index(i)], -coeff / 4))
207 sparse_list.append(('Z', [variables.index(j)], -coeff / 4))
208 new_offset += coeff / 4
209
210 sparse_list.append(('I', [0], new_offset))
211
212 return Hamiltonian(SparsePauliOp.from_sparse_list(sparse_list, num_qubits=len(variables)).simplify())
213
214
[docs]
215class Molecule(Model):
216 def __init__(self, instance: MoleculeInfo, mapper: QubitMapper | None = None, basis_set: str = 'STO-6G') -> None:
217 self.instance = instance
218 self.basis_set = basis_set
219 self.mapper = ParityMapper() if mapper is None else mapper
220 self.problem: ElectronicStructureProblem = PySCFDriver.from_molecule(instance, basis=self.basis_set).run()
221 self.mapper.num_particles = self.problem.num_particles
222 operator = self.mapper.map(self.problem.hamiltonian.second_q_op())
223 if not isinstance(operator, SparsePauliOp):
224 raise TypeError
225 self.operator: SparsePauliOp = operator
226
[docs]
227 @staticmethod
228 def from_preset(instance_name: str) -> 'Molecule':
229 match instance_name:
230 case 'H2':
231 instance = MoleculeInfo(['H', 'H'], [(0.0, 0.0, 0.0), (0.0, 0.0, 0.74)])
232 case _:
233 raise ValueError(f"Molecule {instance_name} not supported, currently you can use: 'H2'")
234 return Molecule(instance)
235
236
[docs]
237class GroverCircuit(Model):
[docs]
238 @staticmethod
239 def create_oracle_from_bitstring(bit_string: str | list[str]) -> Gate:
240 """
241 Creates oracle from given bit string
242 """
243 num_qubits = len(bit_string) + 1
244 qc = QuantumCircuit(num_qubits, name=f'Oracle_{bit_string}')
245
246 reversed_s = bit_string[::-1] # Does this stay???
247
248 qc.x(-1)
249 qc.h(-1)
250
251 for i, char in enumerate(reversed_s):
252 if char == '0':
253 qc.x(i)
254
255 mcx = XGate().control(num_qubits - 1)
256 qc.append(mcx, range(num_qubits))
257
258 for i, char in enumerate(reversed_s):
259 if char == '0':
260 qc.x(i)
261
262 qc.h(-1)
263 qc.x(-1)
264
265 return qc.to_gate()
266
267 @staticmethod
268 def _create_hadamard_walsh_transform(num_qubits: int) -> Gate:
269 print('state_prep is None. Using Hadamard-Walsh Transform')
270 state_prep_circ = QuantumCircuit(num_qubits)
271 state_prep_circ.h(range(num_qubits))
272 return state_prep_circ.to_gate()
273
274 @staticmethod
275 def _validate_and_create_oracle(val: str | list[str]) -> Gate:
276 if not all(c in '01' for c in val):
277 raise ValueError('String/List must contain only zeros and ones.')
278 return GroverCircuit.create_oracle_from_bitstring(val)
279
280 def __init__(
281 self,
282 oracle: QuantumCircuit | np.ndarray | list[str] | str,
283 num_solutions: int = None,
284 num_iterations: int = None,
285 state_prep: QuantumCircuit | Gate | np.ndarray = None,
286 num_state_qubits: int = None,
287 ):
288 """ """
289 self.num_solutions = num_solutions
290
291 if not (num_solutions or num_iterations):
292 raise ValueError('At least one of num_solutions, num_iterations has to be not None')
293
294 self.state_prep = state_prep
295
296 self._gate = None
297
298 oracle_conversion_map = {
299 QuantumCircuit: lambda x: x.to_gate(),
300 str: GroverCircuit._validate_and_create_oracle,
301 list: GroverCircuit._validate_and_create_oracle,
302 np.ndarray: lambda x: UnitaryGate(x),
303 }
304
305 if type(oracle) in oracle_conversion_map:
306 self.oracle = oracle_conversion_map[type(oracle)](oracle)
307 else:
308 raise TypeError(f'Unsupported data type: {type(oracle)}')
309
310 theta = np.arcsin(np.sqrt(num_solutions / (2 ** (self.oracle.num_qubits - 1))))
311 self.num_iterations = num_iterations if num_iterations is not None else int(np.round(np.pi / (4 * theta)))
312
313 self.num_state_qubits = num_state_qubits if num_state_qubits is not None else self.oracle.num_qubits - 1
314
315 state_prep_conversion_map = {
316 QuantumCircuit: lambda x: x.to_gate(),
317 np.ndarray: lambda x: UnitaryGate(x),
318 type(None): lambda _: GroverCircuit._create_hadamard_walsh_transform(self.num_state_qubits),
319 }
320
321 if type(self.state_prep) in state_prep_conversion_map:
322 self.state_prep = state_prep_conversion_map[type(self.state_prep)](self.state_prep)
323
324 self.num_qubits = self.oracle.num_qubits
325
326
[docs]
327class ShorCircuit(Model):
[docs]
328 @staticmethod
329 def create_modular_multiplier_gate_with_uncomputation(factor: int, modulo: int) -> Gate:
330 n_modulo_number_qubits = int(np.floor(np.log2(modulo))) + 1
331
332 def create_adder_gate(factor: int, modulo: int) -> Gate:
333 phase_adder_circuit = QuantumCircuit(n_modulo_number_qubits + 1) # + 1 to avoid overflow
334
335 for qubit_ix in range(n_modulo_number_qubits + 1):
336 phi = 2 * np.pi * factor / 2 ** (n_modulo_number_qubits + 1 - qubit_ix)
337 phase_adder_circuit.rz(phi, qubit_ix)
338
339 return phase_adder_circuit.to_gate()
340
341 def create_modular_adder_gate(factor: int, modulo: int) -> Gate:
342 factor_adder = create_adder_gate(factor, modulo)
343 modulo_adder = create_adder_gate(modulo, modulo)
344
345 controlled_modulo_adder = modulo_adder.control(1)
346
347 factor_substractor = factor_adder.inverse()
348 modulo_substractor = modulo_adder.inverse()
349
350 qft = QFTGate(num_qubits=n_modulo_number_qubits + 1)
351 qft_i = qft.inverse()
352
353 y = QuantumRegister(n_modulo_number_qubits + 1)
354 ancilla = QuantumRegister(1)
355 controls = QuantumRegister(2)
356
357 modular_adder_circuit = QuantumCircuit(controls, y, ancilla)
358
359 modular_adder_circuit.append(factor_adder.control(2), controls[:] + y[:])
360 modular_adder_circuit.append(modulo_substractor, y)
361 modular_adder_circuit.append(qft_i, y)
362 modular_adder_circuit.cx(y[-1], ancilla)
363 modular_adder_circuit.append(qft, y)
364 modular_adder_circuit.append(controlled_modulo_adder, ancilla[:] + y[:])
365 modular_adder_circuit.append(factor_substractor.control(2), controls[:] + y[:])
366 modular_adder_circuit.append(qft_i, y)
367 modular_adder_circuit.x(y[-1])
368 modular_adder_circuit.cx(y[-1], ancilla)
369 modular_adder_circuit.x(y[-1])
370 modular_adder_circuit.append(qft, y[:])
371 modular_adder_circuit.append(factor_adder.control(2), controls[:] + y[:])
372
373 return modular_adder_circuit.to_gate()
374
375 def create_modular_multiplier_gate(factor: int, modulo: int) -> Gate:
376 qft = QFTGate(num_qubits=n_modulo_number_qubits + 1)
377 qft_i = qft.inverse()
378
379 x = QuantumRegister(n_modulo_number_qubits)
380 y = QuantumRegister(n_modulo_number_qubits + 1)
381 ancilla = QuantumRegister(1)
382 control = QuantumRegister(1)
383
384 modular_multiplier_circuit = QuantumCircuit(control, x, y, ancilla)
385
386 modular_multiplier_circuit.append(qft, y)
387
388 for qubit_ix in range(n_modulo_number_qubits):
389 controlled_modular_adder_gate = create_modular_adder_gate(((2**qubit_ix * factor) % modulo), modulo)
390 modular_multiplier_circuit.append(controlled_modular_adder_gate, control[:] + [x[qubit_ix]] + y[:] + ancilla[:])
391
392 modular_multiplier_circuit.append(qft_i, y)
393
394 return modular_multiplier_circuit.to_gate()
395
396 controlled_modular_multiplier_gate = create_modular_multiplier_gate(factor, modulo)
397 controlled_modular_multiplier_gate_i = create_modular_multiplier_gate(pow(factor, -1, modulo), modulo).inverse()
398
399 control = QuantumRegister(1)
400 x = QuantumRegister(n_modulo_number_qubits)
401 y_with_ancilla = QuantumRegister(n_modulo_number_qubits + 1 + 1)
402
403 controlled_modular_multiplier_gate_circuit = QuantumCircuit(control, x, y_with_ancilla)
404
405 controlled_modular_multiplier_gate_circuit.append(controlled_modular_multiplier_gate, control[:] + x[:] + y_with_ancilla[:])
406
407 controlled_swap_gate = SwapGate().control()
408
409 for qubit_ix in range(n_modulo_number_qubits):
410 controlled_modular_multiplier_gate_circuit.append(
411 controlled_swap_gate, control[:] + [x[qubit_ix]] + [y_with_ancilla[qubit_ix]]
412 ) # ancillas are 0 so it's ok
413
414 controlled_modular_multiplier_gate_circuit.append(controlled_modular_multiplier_gate_i, control[:] + x[:] + y_with_ancilla[:])
415
416 return controlled_modular_multiplier_gate_circuit.to_gate()
417
418 def __init__(
419 self,
420 factor: int,
421 modulo: int,
422 n_phase_kickback_qubits: int = None,
423 circuits: list[QuantumCircuit | np.ndarray | Gate] = None,
424 eigen_state_prep: QuantumCircuit | np.ndarray | Gate = None,
425 ):
426 self.factor = factor
427 self.modulo = modulo
428 self.gates = []
429
430 if not n_phase_kickback_qubits:
431 self.n_phase_kickback_qubits = 2 * int(np.floor(np.log2(self.modulo) + 1))
432 else:
433 self.n_phase_kickback_qubits = n_phase_kickback_qubits
434
435 conversion_map = {
436 QuantumCircuit: lambda c: c.to_gate(),
437 np.ndarray: lambda c: UnitaryGate(c),
438 Gate: lambda c: c,
439 type(None): lambda c: None,
440 }
441
442 if circuits is None:
443 for exp in range(2 * (int(np.floor(np.log2(modulo))) + 1)):
444 tmp_factor = (factor ** (2**exp)) % modulo
445 self.gates.append(ShorCircuit.create_modular_multiplier_gate_with_uncomputation(tmp_factor, modulo))
446 else:
447 for circuit in circuits:
448 target_type = type(circuit)
449
450 if target_type in conversion_map:
451 self.gates.append(conversion_map[target_type](circuit))
452 else:
453 raise TypeError(f'Unsupported data type: {target_type}')
454
455 self.num_qubits = self.gates[0].num_qubits
456
457 self.eigen_state_prep = eigen_state_prep
458
459 target_type = type(eigen_state_prep)
460 if target_type in conversion_map:
461 self.eigen_state_prep = conversion_map[target_type](self.eigen_state_prep)
462 else:
463 raise TypeError(f'Unsupported data type: {target_type}')
