Coverage for src/quafel/pipelines/data_generation/nodes.py: 74%
238 statements
« prev ^ index » next coverage.py v7.6.1, created at 2024-09-27 09:52 +0000
1from qiskit.circuit import (
2 ClassicalRegister,
3 QuantumCircuit,
4 CircuitInstruction,
5 ParameterVector,
6 Reset,
7)
8from qiskit.circuit.library import standard_gates
9from qiskit.circuit.exceptions import CircuitError
10from qiskit.quantum_info import partial_trace, random_unitary
11from qiskit import execute
12from qiskit_aer import StatevectorSimulator
13import numpy as np
15import jax
16import jax.numpy as jnp
17from typing import List, Dict, Tuple
18import pandas as pd
19import logging
20import os
21import hashlib
23log = logging.getLogger(__name__)
26def log_circuit(qasm_circuit):
27 return {"circuit_image": None}
30def generate_evaluation_partitions(evaluation_matrix, skip_combinations):
31 partitions = {}
32 idx = 0
33 for f in evaluation_matrix["frameworks"]:
34 if "qubits" in skip_combinations:
35 q = max(evaluation_matrix["qubits"])
36 for d in evaluation_matrix["depths"]:
37 for s in evaluation_matrix["shots"]:
38 partitions[f"{idx}"] = {
39 "framework": f,
40 "qubits": q,
41 "depth": d,
42 "shots": s,
43 }
44 idx += 1
45 elif "depth" in skip_combinations:
46 d = max(evaluation_matrix["depth"])
47 for q in evaluation_matrix["qubits"]:
48 for s in evaluation_matrix["shots"]:
49 partitions[f"{idx}"] = {
50 "framework": f,
51 "qubits": q,
52 "depth": d,
53 "shots": s,
54 }
55 idx += 1
56 elif "shots" in skip_combinations:
57 s = max(evaluation_matrix["shots"])
58 for d in evaluation_matrix["depths"]:
59 for q in evaluation_matrix["qubits"]:
60 partitions[f"{idx}"] = {
61 "framework": f,
62 "qubits": q,
63 "depth": d,
64 "shots": s,
65 }
66 idx += 1
67 else:
68 for q in evaluation_matrix["qubits"]:
69 for d in evaluation_matrix["depths"]:
70 for s in evaluation_matrix["shots"]:
71 partitions[f"{idx}"] = {
72 "framework": f,
73 "qubits": q,
74 "depth": d,
75 "shots": s,
76 }
77 idx += 1
79 # Add row index, because the order of the order
80 # of the rows might get scrambled otherwise
81 eval_partitions = pd.DataFrame(
82 partitions, index=["framework", "qubits", "depth", "shots"]
83 )
84 log.info(f"Generated {eval_partitions.shape[1]} partitions")
85 return {"evaluation_partitions": eval_partitions}
88def extract_partition_data(partition):
89 # TODO: improve this by accessing data by name
90 framework = partition[partition.columns[0]][0]
91 qubits = int(partition[partition.columns[0]][1])
92 depth = int(partition[partition.columns[0]][2])
93 shots = int(partition[partition.columns[0]][3])
95 return {
96 "qubits": qubits,
97 "depth": depth,
98 "n_shots": shots,
99 "framework": framework,
100 }
103def qasm_circuit_to_qiskit(qasm_circuit):
104 qc = QuantumCircuit.from_qasm_str(qasm_circuit)
106 return {
107 "qiskit_circuit": qc,
108 }
111def _random_circuit(
112 num_qubits,
113 depth,
114 max_operands=2,
115 conditional=False,
116 reset=False,
117 seed=None,
118):
119 """
120 <<<<<<< HEAD
121 Code partly taken from Qiskit:
122 qiskit/circuit/random/utils.py
124 Generates a random quantum circuit with the number of qubits
125 and circuit depth provided.
127 The generating method will favor multi-qubit gates.
129 Parameterized gates will have UNBOUND parameters, thus it is
130 required to bind them before the circuit can be evaluated.
132 Args:
133 num_qubits (int): Number of qubits.
134 depth (int): Depth of the circuit.
135 max_operands (int, optional): Maximum number of operands. Defaults to 2.
136 conditional (bool, optional): Cond. Circuit?. Defaults to False.
137 reset (bool, optional): Whether to reset the circuit. Defaults to False.
138 seed (int, optional): Seed for random number generation. Defaults to None.
140 Returns:
141 QuantumCircuit: The generated quantum circuit.
142 """
143 if num_qubits == 0:
144 return QuantumCircuit()
145 if max_operands < 1 or max_operands > 4:
146 raise CircuitError("max_operands must be between 1 and 4")
147 max_operands = max_operands if num_qubits > max_operands else num_qubits
149 gates_1q = [
150 # (Gate class, number of qubits, number of parameters)
151 (standard_gates.IGate, 1, 0),
152 (standard_gates.XGate, 1, 0),
153 (standard_gates.RZGate, 1, 1),
154 (standard_gates.HGate, 1, 0),
155 (standard_gates.RXGate, 1, 1),
156 (standard_gates.RYGate, 1, 1),
157 (standard_gates.SGate, 1, 0),
158 (standard_gates.TGate, 1, 0),
159 (standard_gates.U2Gate, 1, 2),
160 (standard_gates.U3Gate, 1, 3),
161 (standard_gates.YGate, 1, 0),
162 (standard_gates.ZGate, 1, 0),
163 ]
164 if reset:
165 gates_1q.append((Reset, 1, 0))
166 gates_2q = [
167 (standard_gates.CXGate, 2, 0),
168 (standard_gates.CZGate, 2, 0),
169 (standard_gates.SwapGate, 2, 0),
170 ]
171 gates_3q = [
172 (standard_gates.CCXGate, 3, 0),
173 ]
175 # add gates to the overall set, depending on the number of operands
176 gates = gates_1q.copy()
177 if max_operands >= 2:
178 gates.extend(gates_2q)
179 if max_operands >= 3:
180 gates.extend(gates_3q)
181 gates = np.array(
182 gates,
183 dtype=[("class", object), ("num_qubits", np.int16), ("num_params", np.int32)],
184 )
186 # generate a numpy array, that we will use later, to fill "gaps"
187 gates_1q = np.array(gates_1q, dtype=gates.dtype)
189 qc = QuantumCircuit(num_qubits)
191 cr = ClassicalRegister(num_qubits, "c")
192 qc.add_register(cr)
194 rng = np.random.default_rng(seed)
196 qubits = np.array(qc.qubits, dtype=object, copy=True)
198 # Apply arbitrary random operations in layers across all qubits.
199 for layer_number in range(depth):
200 # We generate all the randomness for the layer in one go,
201 # to avoid many separate calls to
202 # the randomisation routines, which can be fairly slow.
204 # This reliably draws too much randomness,
205 # but it's less expensive than looping over more
206 # calls to the rng. After, trim it down by finding the point
207 # when we've used all the qubits.
208 gate_specs = rng.choice(gates, size=len(qubits))
209 cumulative_qubits = np.cumsum(gate_specs["num_qubits"], dtype=np.int16)
210 # Efficiently find the point in the list where the total gates
211 # would use as many as
212 # possible of, but not more than, the number of qubits in the layer.
213 # If there's slack, fill
214 # it with 1q gates.
215 # This will, with favor multi qubit gates, but does not ensure >1 qubit gates
216 # being used, especially with a low number of layers
217 max_index = np.searchsorted(cumulative_qubits, num_qubits, side="right")
218 gate_specs = gate_specs[:max_index]
219 slack = num_qubits - cumulative_qubits[max_index - 1]
220 if slack:
221 gate_specs = np.hstack((gate_specs, rng.choice(gates_1q, size=slack)))
223 # For efficiency in the Python loop, this uses Numpy vectorisation to
224 # pre-calculate the
225 # indices into the lists of qubits and parameters for every gate,
226 # and then suitably
227 # randomises those lists.
228 q_indices = np.empty(len(gate_specs) + 1, dtype=np.int16)
229 p_indices = np.empty(len(gate_specs) + 1, dtype=np.int16)
230 q_indices[0] = p_indices[0] = 0
231 np.cumsum(gate_specs["num_qubits"], out=q_indices[1:])
232 np.cumsum(gate_specs["num_params"], out=p_indices[1:])
233 # parameters = rng.uniform(0, 2 * np.pi, size=p_indices[-1])
234 parameters = ParameterVector(f"p_{layer_number}", p_indices[-1])
235 rng.shuffle(qubits)
237 # We've now generated everything we're going to need.
238 # Now just to add everything. The
239 # conditional check is outside the two loops to make
240 # the more common case of no conditionals
241 # faster, since in Python we don't have a compiler to do this for us.
242 if conditional and layer_number != 0:
243 is_conditional = rng.random(size=len(gate_specs)) < 0.1
244 condition_values = rng.integers(
245 0, 1 << min(num_qubits, 63), size=np.count_nonzero(is_conditional)
246 )
247 c_ptr = 0
248 for gate, q_start, q_end, p_start, p_end, is_cond in zip(
249 gate_specs["class"],
250 q_indices[:-1],
251 q_indices[1:],
252 p_indices[:-1],
253 p_indices[1:],
254 is_conditional,
255 ):
256 operation = gate(*parameters[p_start:p_end])
257 if is_cond:
258 qc.measure(qc.qubits, cr)
259 # The condition values are required to be bigints,
260 # not Numpy's fixed-width type.
261 operation.condition = (cr, int(condition_values[c_ptr]))
262 c_ptr += 1
263 qc._append(
264 CircuitInstruction(
265 operation=operation, qubits=qubits[q_start:q_end]
266 )
267 )
268 else:
269 for gate, q_start, q_end, p_start, p_end in zip(
270 gate_specs["class"],
271 q_indices[:-1],
272 q_indices[1:],
273 p_indices[:-1],
274 p_indices[1:],
275 ):
276 operation = gate(*parameters[p_start:p_end])
277 qc._append(
278 CircuitInstruction(
279 operation=operation, qubits=qubits[q_start:q_end]
280 )
281 )
283 return qc
286def generate_random_qasm_circuit_from_partition(partition, seed):
287 partition_data = extract_partition_data(partition)
288 return {
289 **generate_random_qasm_circuit(
290 partition_data["qubits"], partition_data["depth"], seed
291 ),
292 "n_shots": partition_data["depth"],
293 "framework": partition_data["framework"],
294 }
297def generate_random_qasm_circuit(
298 qubits: int, depth: int, seed: int
299) -> Dict[str, List[float]]:
300 """
301 Generate a random quantum circuit and its representation as a QASM string.
302 Note that the QASM circuit differs from the original circuit in such a way
303 that all qubits are being measured.
305 Args:
306 qubits: Number of qubits.
307 depth: Number of gate layers.
308 seed: Random seed.
310 Returns:
311 A dictionary with the key 'qasm_circuit' containing the QASM string and
312 the key 'circuit' containing the Qiskit circuit object.
313 """
314 qc = _random_circuit(qubits, depth, max_operands=2, seed=seed)
316 rng = np.random.default_rng(seed)
318 # generate all the parameters of the circuit in one go
319 parameter_values = rng.uniform(0, 2 * np.pi, size=len(qc.parameters))
321 # bind the parameters to the circuit
322 bound_circuit = qc.assign_parameters(
323 {p: v for p, v in zip(qc.parameters, parameter_values)}
324 )
325 # measure all of the bound circuit
326 # note that we explicitly NOT use measure_all because this would
327 # add a barrier in the qasm representation that is not recognized by some
328 # frameworks when translating back later
329 bound_circuit.measure(bound_circuit.qubits, *bound_circuit.cregs)
331 # return the bound circuit and the parameterizable circuit
332 return {"qasm_circuit": bound_circuit.qasm(), "qiskit_circuit": qc}
335def generate_evaluation_matrix(
336 min_qubits: int,
337 max_qubits: int,
338 qubits_increment: int,
339 qubits_type: int,
340 min_depth: int,
341 max_depth: int,
342 depth_increment: int,
343 depth_type: int,
344 min_shots: int,
345 max_shots: int,
346 shots_increment: int,
347 shots_type: int,
348 frameworks: List[str],
349):
350 def generate_ticks(min_t, max_t, inc_t, type_t="linear"):
351 if type_t == "linear":
352 ticks = [i for i in range(min_t, max_t + inc_t, inc_t)]
353 elif "exp" in type_t:
354 base = int(type_t.split("exp")[1])
355 ticks = [base**i for i in range(min_t, max_t + inc_t, inc_t)]
356 else:
357 raise ValueError("Unknown base specified and type is not linear")
359 return ticks
361 qubits = generate_ticks(min_qubits, max_qubits, qubits_increment, qubits_type)
362 depths = generate_ticks(min_depth, max_depth, depth_increment, depth_type)
363 shots = generate_ticks(min_shots, max_shots, shots_increment, shots_type)
365 frameworks = frameworks
367 return {
368 "evaluation_matrix": {
369 "frameworks": frameworks,
370 "qubits": qubits,
371 "depths": depths,
372 "shots": shots,
373 }
374 }
377def get_pqc_statevector(
378 circuit: QuantumCircuit, params: np.ndarray, backend, precision: int, cache=True
379) -> float:
380 # Note that this is a jax rng, so it does not matter if we call that multiple times
381 bound_circuit = circuit.bind_parameters(params)
383 # the qasm representation contains the bound parameters, thus it is ok to hash that
384 hs = hashlib.md5(bound_circuit.qasm().encode()).hexdigest()
386 if cache:
387 name = f"pqc_{hs}.npy"
389 cache_folder = ".cache"
390 if not os.path.exists(cache_folder):
391 os.mkdir(cache_folder)
393 file_path = os.path.join(cache_folder, name)
395 if os.path.isfile(file_path):
396 log.debug(f"Loading PQC statevector from {file_path}")
397 loaded_array = np.load(file_path)
398 return loaded_array
400 # execute the PQC circuit with the current set of parameters
401 # ansatz = circuit(params, circuit.num_qubits)
402 result = execute(bound_circuit, backend=backend).result()
404 # extract the statevector from the simulation result
405 U = result.get_statevector(bound_circuit, decimals=precision).data
407 if cache:
408 log.debug(f"Cached PQC statevector into {file_path}")
409 np.save(file_path, U)
411 return U
414def calculate_entangling_capability(
415 circuit: QuantumCircuit,
416 samples_per_parameter: int,
417 samples_per_qubit: int,
418 seed: int,
419) -> Dict[str, float]:
420 """
421 Calculate the entangling capability of a quantum circuit.
422 The strategy is taken from https://doi.org/10.48550/arXiv.1905.10876
423 Implementation inspiration from
424 https://obliviateandsurrender.github.io/blogs/expr.html
426 Sanity Check; The following circuit should yield an entangling capability of 1.
427 qc = QuantumCircuit(2)
428 qc.h(0)
429 qc.cx(0,1)
430 qc.rx(*ParameterVector("x", 1), 0)
432 Note that if qubits are being measured, this will return almost zero
433 because the statevector simulator cannot work on collapsed states.
435 If the circuit doesn't contain parameterizable operations, nan is returned.
437 Args:
438 circuit (QuantumCircuit): The quantum circuit.
439 samples_per_parameter (int): The number of samples to generate.
440 seed (int): The seed for the random number generator.
442 Returns:
443 dict: A dictionary containing the entangling capability value.
444 """
446 def meyer_wallach(circuit, samples, params_shape, precision, rng):
447 if circuit.num_qubits == 1:
448 return 0
450 mw_measure = np.zeros(samples)
452 # FIXME: unify the range for parameters in the circuit
453 # generation method and the sampling here
454 params = rng.uniform(0, 2 * np.pi, (samples, params_shape))
456 # generate a list from [0..num_qubits-1]
457 # we need that later to trace out the corresponding qubits
458 qb = list(range(circuit.num_qubits))
460 backend = StatevectorSimulator(precision="single")
461 backend.set_options(
462 max_parallel_threads=10,
463 max_parallel_experiments=0,
464 statevector_parallel_threshold=16,
465 )
467 # outer sum of the MW measure; iterate over set of parameters
468 for i in range(samples):
469 U = get_pqc_statevector(circuit, params[i], backend, precision)
471 # generate a list from [0..num_qubits-1]
472 # we need that later to trace out the corresponding qubits
473 qb = list(range(circuit.num_qubits))
474 # initialize the inner sum which corresponds to the entropy
475 entropy = 0
477 # inner sum of the MW measure
478 for j in range(circuit.num_qubits):
479 # density of the jth qubit after tracing out the rest
480 density = partial_trace(U, qb[:j] + qb[j + 1 :]).data
481 # trace of the density matrix
482 entropy += np.trace(density @ density).real
484 # fixes accumulating decimals that would otherwise lead to a MW > 1
485 entropy = entropy / circuit.num_qubits
486 # inverse of the normalized entropy is the MW
487 # for the current sample of parameters
488 mw_measure[i] = 1 - entropy
489 # final normalization according to formula
490 return max(0, min((2 * np.sum(mw_measure) / samples), 1))
492 rng = np.random.default_rng(seed=seed)
494 if len(circuit.parameters) == 0:
495 entangling_capability = 0
496 else:
497 # TODO: propagate precision to kedro parameters
498 entangling_capability = meyer_wallach(
499 circuit=circuit,
500 samples=samples_per_parameter * int(np.log(len(circuit.parameters)) + 1)
501 + samples_per_qubit * circuit.num_qubits,
502 params_shape=len(circuit.parameters),
503 precision=5,
504 rng=rng,
505 )
507 return {"entangling_capability": entangling_capability}
510def calculate_expressibility(
511 circuit: QuantumCircuit,
512 samples_per_parameter: int,
513 haar_samples_per_qubit: int,
514 seed: int,
515) -> Dict[str, float]:
516 """
517 Calculate the expressibility of a PQC circuit using
518 a randomized estimation scheme.
519 The strategy is taken from https://doi.org/10.48550/arXiv.1905.10876
520 Implementation inspiration from
521 https://obliviateandsurrender.github.io/blogs/expr.html
523 Args:
524 circuit (QuantumCircuit): The PQC circuit to be analyzed
525 samples_per_parameter (int): The number of samples to use for estimation
526 seed (int): The seed for the random number generator
528 Returns:
529 Dict[str, float]: A dictionary containing the expressibility value
530 """
532 def random_haar_unitary(n_qubits: int, rng: np.random.RandomState) -> np.ndarray:
533 """
534 Generate a random unitary matrix in the Haar measure.
535 For details on the QR decomposition, see
536 https://doi.org/10.48550/arXiv.math-ph/0609050
538 Args:
539 n_qubits (int): The number of qubits in the system
540 rng (np.random.RandomState): The RandomState object to use for
541 generating random numbers
543 Returns:
544 np.ndarray: A 2^n x 2^n unitary matrix representing a random
545 unitary in the Haar measure on the unitary group U(2^n)
546 """
548 N = 2**n_qubits
550 # # Generate uniformly sampled random complex numbers
551 # Z = jax.random.normal(rng, (N, N)) + 1.0j * jax.random.normal(rng, (N, N))
553 # def f(Z):
554 # # Do a QR decomposition
555 # [Q, R] = jnp.linalg.qr(Z)
556 # # Allow D following unitary matrix constraints
557 # D = jnp.diag(jnp.diagonal(R) / jnp.abs(jnp.diagonal(R)))
558 # # Composite the Haar Unitary
559 # return jnp.dot(Q, D)
560 # return jax.jit(f)(Z)
562 return random_unitary(
563 N,
564 int(rng._base_array.real[0]),
565 ).data
567 def haar_integral(
568 n_qubits: int, samples: int, rng: np.random.RandomState
569 ) -> np.ndarray:
570 """
571 Compute the Haar integral for a given number of samples
573 Args:
574 n_qubits (int): The number of qubits in the system
575 samples (int): The number of samples to use for estimation
576 rng (np.random.RandomState): The RandomState object to use for
577 generating random numbers
579 Returns:
580 np.ndarray: A 2^n x 2^n array representing the normalized max
581 expressibility
582 """
583 N = 2**n_qubits
584 Z = jnp.zeros((N, N), dtype=complex)
586 zero_state = jnp.zeros(N, dtype=complex)
587 zero_state = zero_state.at[0].set(1)
589 def f(zero_state, U):
590 A = jnp.matmul(zero_state, U).reshape(-1, 1)
591 return jnp.kron(A, A.conj().T)
593 jf = jax.jit(f)
595 for i in range(samples):
596 rng, subkey = jax.random.split(rng)
597 U = random_haar_unitary(n_qubits, subkey)
599 Z += jf(zero_state, U)
601 return Z / samples
603 def get_haar_integral(
604 n_qubits: int, samples: int, rng: np.random.RandomState, cache=True
605 ) -> float:
606 if cache:
607 name = f"haar_{n_qubits}q_{samples}s.npy"
609 cache_folder = ".cache"
610 if not os.path.exists(cache_folder):
611 os.mkdir(cache_folder)
613 file_path = os.path.join(cache_folder, name)
615 if os.path.isfile(file_path):
616 log.debug(f"Loading Haar integral from {file_path}")
617 loaded_array = np.load(file_path)
618 return loaded_array
620 # Note that this is a jax rng, so it does not matter if we
621 # call that multiple times
622 array = haar_integral(n_qubits, samples, rng)
624 if cache:
625 log.debug(f"Cached Haar integral into {file_path}")
626 np.save(file_path, array)
628 return array
630 def pqc_integral(
631 circuit: QuantumCircuit,
632 samples: int,
633 params_shape: Tuple,
634 precision: int,
635 rng: np.random.default_rng,
636 ) -> np.ndarray:
637 """
638 Compute the entangling capability of a PQC circuit using a randomized
639 estimation scheme
641 Args:
642 circuit (QuantumCircuit): The PQC circuit to be analyzed
643 samples (int): The number of samples to use for estimation
644 params_shape (tuple): The shape of the array of parameters to be
645 used in the circuit simulation
646 precision (int): The number of decimals to use when extracting the
647 statevector from the simulation result
649 Returns:
650 np.ndarray: A 2^n x 2^n array representing the expressibility
651 of the PQC circuit, where n is the number of qubits in the circuit
652 """
653 N = 2**circuit.num_qubits
654 Z = jnp.zeros((N, N), dtype=complex)
656 def f(U):
657 return jnp.kron(U, U.conj().T) # type: ignore
659 jf = jax.jit(f)
661 backend = StatevectorSimulator(precision="single")
662 backend.set_options(
663 max_parallel_threads=10,
664 max_parallel_experiments=0,
665 statevector_parallel_threshold=16,
666 )
668 # FIXME: unify the range for parameters in the circuit
669 # generation method and the sampling here
670 params = rng.uniform(0, 2 * np.pi, (samples, params_shape))
671 for i in range(samples):
672 # extract the statevector from the simulation result
673 U = get_pqc_statevector(circuit, params[i], backend, precision).reshape(
674 -1, 1
675 )
677 # accumulate the contribution to the expressibility
678 Z += jf(U)
680 return Z / samples
682 rng = np.random.default_rng(seed=seed)
683 jrng = jax.random.key(seed)
685 def f(haar, pqc):
686 # Note that we use the INVERSE here, because a
687 # a LOW distance would actually correspond to a HIGH expressibility
688 return 1 - jnp.linalg.norm(haar - pqc)
690 jf = jax.jit(f)
692 # FIXME: the actual value is strongly dependend on the seed (~5-10% deviation)
693 # TODO: propagate precision to kedro parameters
694 if len(circuit.parameters) == 0:
695 expressibility = 0
696 else:
697 hi = get_haar_integral(
698 n_qubits=circuit.num_qubits, samples=haar_samples_per_qubit, rng=jrng
699 )
700 pi = pqc_integral(
701 circuit=circuit,
702 samples=samples_per_parameter * int(np.log(len(circuit.parameters)) + 1)
703 + haar_samples_per_qubit * circuit.num_qubits,
704 params_shape=len(circuit.parameters),
705 precision=5,
706 rng=rng,
707 )
708 expressibility = jf(hi, pi)
710 return {
711 "expressibility": expressibility,
712 }
715def calculate_measures(
716 circuit: QuantumCircuit,
717 samples_per_parameter: int,
718 haar_samples_per_qubit: int,
719 seed: int,
720) -> List[float]:
721 expressibility = calculate_expressibility(
722 circuit=circuit,
723 samples_per_parameter=samples_per_parameter,
724 haar_samples_per_qubit=haar_samples_per_qubit,
725 seed=seed,
726 )["expressibility"]
728 entangling_capability = calculate_entangling_capability(
729 circuit=circuit,
730 samples_per_parameter=samples_per_parameter,
731 samples_per_qubit=haar_samples_per_qubit,
732 seed=seed,
733 )["entangling_capability"]
735 return {
736 "measure": pd.DataFrame(
737 {
738 "expressibility": [expressibility],
739 "entangling_capability": [entangling_capability],
740 }
741 )
742 }