6835383038

Committed 11 Nov 2023 04:33PM UTC coverage: 85.914% (-1.0%) from 86.937%

Build # 6835383038

Build Type

push

github

Committed by

web-flow

Commit Message

Remove deprecated ``qiskit.algorithms`` module (#11086)

* Delete qiskit.algorithms and tests

* Remove algorithms from tests

* Remove references from docs

* Remove references from renos

* Remove optimizer references from AQC

* Remove algorithms API doc

* Add release note

* Revert "Remove references from renos"

This reverts commit 6b291ab04.

* Update releasenotes/notes/remove-qiskit-algorithms-a43541fe24b72208.yaml

Run Details

65911 of 76717 relevant lines covered (85.91%)

100402.81 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

96.97

/qiskit/quantum_info/operators/symplectic/random.py

# This code is part of Qiskit.
#
# (C) Copyright IBM 2017, 2020
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.
"""
Random symplectic operator functions
"""

from __future__ import annotations

import numpy as np
from numpy.random import default_rng

from .clifford import Clifford
from .pauli import Pauli
from .pauli_list import PauliList


def random_pauli(
    num_qubits: int, group_phase: bool = False, seed: int | np.random.Generator | None = None
):
    """Return a random Pauli.

    Args:
        num_qubits (int): the number of qubits.
        group_phase (bool): Optional. If True generate random phase.
                            Otherwise the phase will be set so that the
                            Pauli coefficient is +1 (default: False).
        seed (int or np.random.Generator): Optional. Set a fixed seed or
                                           generator for RNG.

    Returns:
        Pauli: a random Pauli
    """
    if seed is None:
        rng = np.random.default_rng()
    elif isinstance(seed, np.random.Generator):
        rng = seed
    else:
        rng = default_rng(seed)
    z = rng.integers(2, size=num_qubits, dtype=bool)
    x = rng.integers(2, size=num_qubits, dtype=bool)
    phase = rng.integers(4) if group_phase else 0
    pauli = Pauli((z, x, phase))
    return pauli


def random_pauli_list(
    num_qubits: int,
    size: int = 1,
    seed: int | np.random.Generator | None = None,
    phase: bool = True,
):
    """Return a random PauliList.

    Args:
        num_qubits (int): the number of qubits.
        size (int): Optional. The length of the Pauli list (Default: 1).
        seed (int or np.random.Generator): Optional. Set a fixed seed or generator for RNG.
        phase (bool): If True the Pauli phases are randomized, otherwise the phases are fixed to 0.
                     [Default: True]

    Returns:
        PauliList: a random PauliList.
    """
    if seed is None:
        rng = np.random.default_rng()
    elif isinstance(seed, np.random.Generator):
        rng = seed
    else:
        rng = default_rng(seed)

    z = rng.integers(2, size=(size, num_qubits)).astype(bool)
    x = rng.integers(2, size=(size, num_qubits)).astype(bool)
    if phase:
        _phase = rng.integers(4, size=(size))
        return PauliList.from_symplectic(z, x, _phase)
    return PauliList.from_symplectic(z, x)


def random_clifford(num_qubits: int, seed: int | np.random.Generator | None = None):
    """Return a random Clifford operator.

    The Clifford is sampled using the method of Reference [1].

    Args:
        num_qubits (int): the number of qubits for the Clifford
        seed (int or np.random.Generator): Optional. Set a fixed seed or
                                           generator for RNG.

    Returns:
        Clifford: a random Clifford operator.

    Reference:
        1. S. Bravyi and D. Maslov, *Hadamard-free circuits expose the
           structure of the Clifford group*.
           `arXiv:2003.09412 [quant-ph] <https://arxiv.org/abs/2003.09412>`_
    """
    if seed is None:
        rng = np.random.default_rng()
    elif isinstance(seed, np.random.Generator):
        rng = seed
    else:
        rng = default_rng(seed)

    had, perm = _sample_qmallows(num_qubits, rng)

    gamma1 = np.diag(rng.integers(2, size=num_qubits, dtype=np.int8))
    gamma2 = np.diag(rng.integers(2, size=num_qubits, dtype=np.int8))
    delta1 = np.eye(num_qubits, dtype=np.int8)
    delta2 = delta1.copy()

    _fill_tril(gamma1, rng, symmetric=True)
    _fill_tril(gamma2, rng, symmetric=True)
    _fill_tril(delta1, rng)
    _fill_tril(delta2, rng)

    # For large num_qubits numpy.inv function called below can
    # return invalid output leading to a non-symplectic Clifford
    # being generated. This can be prevented by manually forcing
    # block inversion of the matrix.
    block_inverse_threshold = 50

    # Compute stabilizer table
    zero = np.zeros((num_qubits, num_qubits), dtype=np.int8)
    prod1 = np.matmul(gamma1, delta1) % 2
    prod2 = np.matmul(gamma2, delta2) % 2
    inv1 = _inverse_tril(delta1, block_inverse_threshold).transpose()
    inv2 = _inverse_tril(delta2, block_inverse_threshold).transpose()
    table1 = np.block([[delta1, zero], [prod1, inv1]])
    table2 = np.block([[delta2, zero], [prod2, inv2]])

    # Apply qubit permutation
    table = table2[np.concatenate([perm, num_qubits + perm])]

    # Apply layer of Hadamards
    inds = had * np.arange(1, num_qubits + 1)
    inds = inds[inds > 0] - 1
    lhs_inds = np.concatenate([inds, inds + num_qubits])
    rhs_inds = np.concatenate([inds + num_qubits, inds])
    table[lhs_inds, :] = table[rhs_inds, :]

    # Apply table
    tableau = np.zeros((2 * num_qubits, 2 * num_qubits + 1), dtype=bool)
    tableau[:, :-1] = np.mod(np.matmul(table1, table), 2)

    # Generate random phases
    tableau[:, -1] = rng.integers(2, size=2 * num_qubits)
    return Clifford(tableau, validate=False)


def _sample_qmallows(n, rng=None):
    """Sample from the quantum Mallows distribution"""

    if rng is None:
        rng = np.random.default_rng()

    # Hadmard layer
    had = np.zeros(n, dtype=bool)

    # Permutation layer
    perm = np.zeros(n, dtype=int)

    inds = list(range(n))
    for i in range(n):
        m = n - i
        eps = 4 ** (-m)
        r = rng.uniform(0, 1)
        index = -int(np.ceil(np.log2(r + (1 - r) * eps)))
        had[i] = index < m
        if index < m:
            k = index
        else:
            k = 2 * m - index - 1
        perm[i] = inds[k]
        del inds[k]
    return had, perm


def _fill_tril(mat, rng, symmetric=False):
    """Add symmetric random ints to off diagonals"""
    dim = mat.shape[0]
    # Optimized for low dimensions
    if dim == 1:
        return

    if dim <= 4:
        mat[1, 0] = rng.integers(2, dtype=np.int8)
        if symmetric:
            mat[0, 1] = mat[1, 0]
        if dim > 2:
            mat[2, 0] = rng.integers(2, dtype=np.int8)
            mat[2, 1] = rng.integers(2, dtype=np.int8)
            if symmetric:
                mat[0, 2] = mat[2, 0]
                mat[1, 2] = mat[2, 1]
        if dim > 3:
            mat[3, 0] = rng.integers(2, dtype=np.int8)
            mat[3, 1] = rng.integers(2, dtype=np.int8)
            mat[3, 2] = rng.integers(2, dtype=np.int8)
            if symmetric:
                mat[0, 3] = mat[3, 0]
                mat[1, 3] = mat[3, 1]
                mat[2, 3] = mat[3, 2]
        return

    # Use numpy indices for larger dimensions
    rows, cols = np.tril_indices(dim, -1)
    vals = rng.integers(2, size=rows.size, dtype=np.int8)
    mat[(rows, cols)] = vals
    if symmetric:
        mat[(cols, rows)] = vals


def _inverse_tril(mat, block_inverse_threshold):
    """Invert a lower-triangular matrix with unit diagonal."""
    # Optimized inversion function for low dimensions
    dim = mat.shape[0]

    if dim <= 2:
        return mat

    if dim <= 5:
        inv = mat.copy()
        inv[2, 0] = mat[2, 0] ^ (mat[1, 0] & mat[2, 1])
        if dim > 3:
            inv[3, 1] = mat[3, 1] ^ (mat[2, 1] & mat[3, 2])
            inv[3, 0] = mat[3, 0] ^ (mat[3, 2] & mat[2, 0]) ^ (mat[1, 0] & inv[3, 1])
        if dim > 4:
            inv[4, 2] = (mat[4, 2] ^ (mat[3, 2] & mat[4, 3])) & 1
            inv[4, 1] = mat[4, 1] ^ (mat[4, 3] & mat[3, 1]) ^ (mat[2, 1] & inv[4, 2])
            inv[4, 0] = (
                mat[4, 0]
                ^ (mat[1, 0] & inv[4, 1])
                ^ (mat[2, 0] & inv[4, 2])
                ^ (mat[3, 0] & mat[4, 3])
            )
        return inv % 2

    # For higher dimensions we use Numpy's inverse function
    # however this function tends to fail and result in a non-symplectic
    # final matrix if n is too large.
    if dim <= block_inverse_threshold:
        return np.linalg.inv(mat).astype(np.int8) % 2

    # For very large matrices  we divide the matrix into 4 blocks of
    # roughly equal size and use the analytic formula for the inverse
    # of a block lower-triangular matrix:
    # inv([[A, 0],[C, D]]) = [[inv(A), 0], [inv(D).C.inv(A), inv(D)]]
    # call the inverse function recursively to compute inv(A) and invD

    dim1 = dim // 2
    mat_a = _inverse_tril(mat[0:dim1, 0:dim1], block_inverse_threshold)
    mat_d = _inverse_tril(mat[dim1:dim, dim1:dim], block_inverse_threshold)
    mat_c = np.matmul(np.matmul(mat_d, mat[dim1:dim, 0:dim1]), mat_a)
    inv = np.block([[mat_a, np.zeros((dim1, dim - dim1), dtype=int)], [mat_c, mat_d]])
    return inv % 2

1	# This code is part of Qiskit.
2	#
3	# (C) Copyright IBM 2017, 2020
4	#
5	# This code is licensed under the Apache License, Version 2.0. You may
6	# obtain a copy of this license in the LICENSE.txt file in the root directory
7	# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
8	#
9	# Any modifications or derivative works of this code must retain this
10	# copyright notice, and modified files need to carry a notice indicating
11	# that they have been altered from the originals.
12	"""	1✔
13	Random symplectic operator functions
14	"""
15
16	from __future__ import annotations	1✔
17
18	import numpy as np	1✔
19	from numpy.random import default_rng	1✔
20
21	from .clifford import Clifford	1✔
22	from .pauli import Pauli	1✔
23	from .pauli_list import PauliList	1✔
24
25
26	def random_pauli(	1✔
27	num_qubits: int, group_phase: bool = False, seed: int \| np.random.Generator \| None = None
28	):
29	"""Return a random Pauli.
30
31	Args:
32	num_qubits (int): the number of qubits.
33	group_phase (bool): Optional. If True generate random phase.
34	Otherwise the phase will be set so that the
35	Pauli coefficient is +1 (default: False).
36	seed (int or np.random.Generator): Optional. Set a fixed seed or
37	generator for RNG.
38
39	Returns:
40	Pauli: a random Pauli
41	"""
42	if seed is None:	1✔
43	rng = np.random.default_rng()	×
44	elif isinstance(seed, np.random.Generator):	1✔
45	rng = seed	1✔
46	else:
47	rng = default_rng(seed)	1✔
48	z = rng.integers(2, size=num_qubits, dtype=bool)	1✔
49	x = rng.integers(2, size=num_qubits, dtype=bool)	1✔
50	phase = rng.integers(4) if group_phase else 0	1✔
51	pauli = Pauli((z, x, phase))	1✔
52	return pauli	1✔
53
54
55	def random_pauli_list(	1✔
56	num_qubits: int,
57	size: int = 1,
58	seed: int \| np.random.Generator \| None = None,
59	phase: bool = True,
60	):
61	"""Return a random PauliList.
62
63	Args:
64	num_qubits (int): the number of qubits.
65	size (int): Optional. The length of the Pauli list (Default: 1).
66	seed (int or np.random.Generator): Optional. Set a fixed seed or generator for RNG.
67	phase (bool): If True the Pauli phases are randomized, otherwise the phases are fixed to 0.
68	[Default: True]
69
70	Returns:
71	PauliList: a random PauliList.
72	"""
73	if seed is None:	1✔
74	rng = np.random.default_rng()	1✔
75	elif isinstance(seed, np.random.Generator):	1✔
76	rng = seed	×
77	else:
78	rng = default_rng(seed)	1✔
79
80	z = rng.integers(2, size=(size, num_qubits)).astype(bool)	1✔
81	x = rng.integers(2, size=(size, num_qubits)).astype(bool)	1✔
82	if phase:	1✔
83	_phase = rng.integers(4, size=(size))	1✔
84	return PauliList.from_symplectic(z, x, _phase)	1✔
85	return PauliList.from_symplectic(z, x)	1✔
86
87
88	def random_clifford(num_qubits: int, seed: int \| np.random.Generator \| None = None):	1✔
89	"""Return a random Clifford operator.
90
91	The Clifford is sampled using the method of Reference [1].
92
93	Args:
94	num_qubits (int): the number of qubits for the Clifford
95	seed (int or np.random.Generator): Optional. Set a fixed seed or
96	generator for RNG.
97
98	Returns:
99	Clifford: a random Clifford operator.
100
101	Reference:
102	1. S. Bravyi and D. Maslov, *Hadamard-free circuits expose the
103	structure of the Clifford group*.
104	`arXiv:2003.09412 [quant-ph] <https://arxiv.org/abs/2003.09412>`_
105	"""
106	if seed is None:	1✔
107	rng = np.random.default_rng()	×
108	elif isinstance(seed, np.random.Generator):	1✔
109	rng = seed	1✔
110	else:
111	rng = default_rng(seed)	1✔
112
113	had, perm = _sample_qmallows(num_qubits, rng)	1✔
114
115	gamma1 = np.diag(rng.integers(2, size=num_qubits, dtype=np.int8))	1✔
116	gamma2 = np.diag(rng.integers(2, size=num_qubits, dtype=np.int8))	1✔
117	delta1 = np.eye(num_qubits, dtype=np.int8)	1✔
118	delta2 = delta1.copy()	1✔
119
120	_fill_tril(gamma1, rng, symmetric=True)	1✔
121	_fill_tril(gamma2, rng, symmetric=True)	1✔
122	_fill_tril(delta1, rng)	1✔
123	_fill_tril(delta2, rng)	1✔
124
125	# For large num_qubits numpy.inv function called below can
126	# return invalid output leading to a non-symplectic Clifford
127	# being generated. This can be prevented by manually forcing
128	# block inversion of the matrix.
129	block_inverse_threshold = 50	1✔
130
131	# Compute stabilizer table
132	zero = np.zeros((num_qubits, num_qubits), dtype=np.int8)	1✔
133	prod1 = np.matmul(gamma1, delta1) % 2	1✔
134	prod2 = np.matmul(gamma2, delta2) % 2	1✔
135	inv1 = _inverse_tril(delta1, block_inverse_threshold).transpose()	1✔
136	inv2 = _inverse_tril(delta2, block_inverse_threshold).transpose()	1✔
137	table1 = np.block([[delta1, zero], [prod1, inv1]])	1✔
138	table2 = np.block([[delta2, zero], [prod2, inv2]])	1✔
139
140	# Apply qubit permutation
141	table = table2[np.concatenate([perm, num_qubits + perm])]	1✔
142
143	# Apply layer of Hadamards
144	inds = had * np.arange(1, num_qubits + 1)	1✔
145	inds = inds[inds > 0] - 1	1✔
146	lhs_inds = np.concatenate([inds, inds + num_qubits])	1✔
147	rhs_inds = np.concatenate([inds + num_qubits, inds])	1✔
148	table[lhs_inds, :] = table[rhs_inds, :]	1✔
149
150	# Apply table
151	tableau = np.zeros((2 * num_qubits, 2 * num_qubits + 1), dtype=bool)	1✔
152	tableau[:, :-1] = np.mod(np.matmul(table1, table), 2)	1✔
153
154	# Generate random phases
155	tableau[:, -1] = rng.integers(2, size=2 * num_qubits)	1✔
156	return Clifford(tableau, validate=False)	1✔
157
158
159	def _sample_qmallows(n, rng=None):	1✔
160	"""Sample from the quantum Mallows distribution"""
161
162	if rng is None:	1✔
163	rng = np.random.default_rng()	×
164
165	# Hadmard layer
166	had = np.zeros(n, dtype=bool)	1✔
167
168	# Permutation layer
169	perm = np.zeros(n, dtype=int)	1✔
170
171	inds = list(range(n))	1✔
172	for i in range(n):	1✔
173	m = n - i	1✔
174	eps = 4 ** (-m)	1✔
175	r = rng.uniform(0, 1)	1✔
176	index = -int(np.ceil(np.log2(r + (1 - r) * eps)))	1✔
177	had[i] = index < m	1✔
178	if index < m:	1✔
179	k = index	1✔
180	else:
181	k = 2 * m - index - 1	1✔
182	perm[i] = inds[k]	1✔
183	del inds[k]	1✔
184	return had, perm	1✔
185
186
187	def _fill_tril(mat, rng, symmetric=False):	1✔
188	"""Add symmetric random ints to off diagonals"""
189	dim = mat.shape[0]	1✔
190	# Optimized for low dimensions
191	if dim == 1:	1✔
192	return	1✔
193
194	if dim <= 4:	1✔
195	mat[1, 0] = rng.integers(2, dtype=np.int8)	1✔
196	if symmetric:	1✔
197	mat[0, 1] = mat[1, 0]	1✔
198	if dim > 2:	1✔
199	mat[2, 0] = rng.integers(2, dtype=np.int8)	1✔
200	mat[2, 1] = rng.integers(2, dtype=np.int8)	1✔
201	if symmetric:	1✔
202	mat[0, 2] = mat[2, 0]	1✔
203	mat[1, 2] = mat[2, 1]	1✔
204	if dim > 3:	1✔
205	mat[3, 0] = rng.integers(2, dtype=np.int8)	1✔
206	mat[3, 1] = rng.integers(2, dtype=np.int8)	1✔
207	mat[3, 2] = rng.integers(2, dtype=np.int8)	1✔
208	if symmetric:	1✔
209	mat[0, 3] = mat[3, 0]	1✔
210	mat[1, 3] = mat[3, 1]	1✔
211	mat[2, 3] = mat[3, 2]	1✔
212	return	1✔
213
214	# Use numpy indices for larger dimensions
215	rows, cols = np.tril_indices(dim, -1)	1✔
216	vals = rng.integers(2, size=rows.size, dtype=np.int8)	1✔
217	mat[(rows, cols)] = vals	1✔
218	if symmetric:	1✔
219	mat[(cols, rows)] = vals	1✔
220
221
222	def _inverse_tril(mat, block_inverse_threshold):	1✔
223	"""Invert a lower-triangular matrix with unit diagonal."""
224	# Optimized inversion function for low dimensions
225	dim = mat.shape[0]	1✔
226
227	if dim <= 2:	1✔
228	return mat	1✔
229
230	if dim <= 5:	1✔
231	inv = mat.copy()	1✔
232	inv[2, 0] = mat[2, 0] ^ (mat[1, 0] & mat[2, 1])	1✔
233	if dim > 3:	1✔
234	inv[3, 1] = mat[3, 1] ^ (mat[2, 1] & mat[3, 2])	1✔
235	inv[3, 0] = mat[3, 0] ^ (mat[3, 2] & mat[2, 0]) ^ (mat[1, 0] & inv[3, 1])	1✔
236	if dim > 4:	1✔
237	inv[4, 2] = (mat[4, 2] ^ (mat[3, 2] & mat[4, 3])) & 1	1✔
238	inv[4, 1] = mat[4, 1] ^ (mat[4, 3] & mat[3, 1]) ^ (mat[2, 1] & inv[4, 2])	1✔
239	inv[4, 0] = (	1✔
240	mat[4, 0]
241	^ (mat[1, 0] & inv[4, 1])
242	^ (mat[2, 0] & inv[4, 2])
243	^ (mat[3, 0] & mat[4, 3])
244	)
245	return inv % 2	1✔
246
247	# For higher dimensions we use Numpy's inverse function
248	# however this function tends to fail and result in a non-symplectic
249	# final matrix if n is too large.
250	if dim <= block_inverse_threshold:	1✔
251	return np.linalg.inv(mat).astype(np.int8) % 2	1✔
252
253	# For very large matrices we divide the matrix into 4 blocks of
254	# roughly equal size and use the analytic formula for the inverse
255	# of a block lower-triangular matrix:
256	# inv([[A, 0],[C, D]]) = [[inv(A), 0], [inv(D).C.inv(A), inv(D)]]
257	# call the inverse function recursively to compute inv(A) and invD
258
259	dim1 = dim // 2	1✔
260	mat_a = _inverse_tril(mat[0:dim1, 0:dim1], block_inverse_threshold)	1✔
261	mat_d = _inverse_tril(mat[dim1:dim, dim1:dim], block_inverse_threshold)	1✔
262	mat_c = np.matmul(np.matmul(mat_d, mat[dim1:dim, 0:dim1]), mat_a)	1✔
263	inv = np.block([[mat_a, np.zeros((dim1, dim - dim1), dtype=int)], [mat_c, mat_d]])	1✔
264	return inv % 2	1✔

Qiskit / qiskit / 6835383038

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous