6893865468

Committed 16 Nov 2023 04:49PM UTC coverage: 43.147% (+0.7%) from 42.467%

Build # 6893865468

Build Type

push

github

Committed by

feihoo87

Commit Message

update

Run Details

7 of 17 new or added lines in 4 files covered. (41.18%)

588 existing lines in 10 files now uncovered.

7436 of 17234 relevant lines covered (43.15%)

3.87 hits per line

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

37.86

/waveforms/quantum/math.py

from functools import reduce
from itertools import chain, product

import numpy as np
from scipy.linalg import eigh, expm, logm, sqrtm

from waveforms.math.matricies import (BellPhiM, BellPhiP, BellPsiM, BellPsiP,
                                      sigmaI, sigmaX, sigmaY, sigmaZ)
from waveforms.math.signal import decay, oscillation

# Paulis
s0, s1, s2, s3 = sigmaI(), sigmaX(), sigmaY(), sigmaZ()

# Bell states
phiplus, phiminus = BellPhiP(), BellPhiM()
psiplus, psiminus = BellPsiP(), BellPsiM()


def tensor(matrixes: np.ndarray) -> np.ndarray:
    """Compute the tensor product of a list (or array) of matrixes"""
    return reduce(np.kron, matrixes)


def sigma(j: int, N: int = 1) -> np.ndarray:
    """
    """
    s = [s0, s1, s2, s3]
    dims = [4] * N
    idx = np.unravel_index(j, dims)
    return tensor(s[x] for x in idx)


def basis(x: int | str, dim: int = 2) -> np.ndarray:
    '''
    Returns a single element of a state vector.  Each component is
    either a digit (0,1,2,...) representing the computational basis,
    or one or two letters representing a special state (H,V,D,A,R,L,
    X,Y,Z,Xm,Ym,Zm.  B1..4 are the 4 bell states
    '''
    d = {
        'H': [1, 0],
        'V': [0, 1],
        'D': [1, 1],
        'A': [1, -1],
        'R': [1, 1j],
        'L': [1, -1j],
        'X': [1, 1],
        'Xp': [1, 1],
        'Xm': [1, -1],
        'Y': [1, 1j],
        'Yp': [1, 1j],
        'Ym': [1, -1j],
        'Z': [1, 0],
        'Zp': [1, 0],
        'Zm': [0, 1],
        'B1': phiminus,
        'B2': phiplus,
        'B3': psiminus,
        'B4': psiplus
    }
    if isinstance(x, int) or x.isdigit():
        rv = np.zeros((dim), dtype=complex)
        rv[int(x)] = 1
    else:
        rv = np.array(d[x], dtype=complex)
        rv = rv / np.linalg.norm(rv)
    return rv


def commutator(A: np.ndarray, B: np.ndarray) -> np.ndarray:
    """Compute the commutator of two matrices"""
    return A @ B - B @ A


def dagger(vector: np.ndarray) -> np.ndarray:
    """Compute the hermitian conjugate of a vector or matrix"""
    return vector.conjugate().transpose()


def normalize(X: np.ndarray) -> np.ndarray:
    """Normalize the representation of X.  This means:
    For state vectors, X' * X = 1 and the first non-zero
    element of X is positive real.  For density matricies,
    X is hermetian, unit trace, and non-negative"""
    if X.ndim == 1:
        nz = X[np.nonzero(X)][0]
        X *= nz.conj()
        amp = np.abs(dagger(X) @ X)
        X /= np.sqrt(amp)
    else:
        (d, U) = eigh(X)
        d = d.real
        d[d < 0] = 0
        X = U @ np.diag(d) @ dagger(U)
        X = dagger(X) + X
        X /= np.trace(X)
    return X


def psi2rho(psi: np.ndarray) -> np.ndarray:
    '''
    Converts input state vector to density matrix represenation.  If psi
    is already a density matrix, leave it alone, so this function can be
    used to coerce input arguments to the desired type.
    '''
    if psi.ndim == 1:
        return normalize(np.outer(psi, psi.conj()))
    else:  # Already a density matrix
        return normalize(psi)


def rho2v(rho):
    """密度矩阵转相干矢
    """
    N = rho.shape[0]
    indices = np.tril_indices(N, -1)
    z = np.diag(rho).real
    Z = np.cumsum(z[:-1]) - z[1:] * np.arange(1, N)
    return np.hstack((rho[indices].real, rho[indices].imag, Z))


def v2rho(V):
    """相干矢转密度矩阵"""
    N = int(np.sqrt(len(V) + 1))
    assert N**2 - 1 == len(V)

    X = V[:(N**2 - N) // 2]
    Y = V[(N**2 - N) // 2:(N**2 - N)]
    Z = V[(N**2 - N):]

    A = np.hstack([Z[:0:-1] - Z[-2::-1], Z[0]])
    zn = (1 - A.sum()) / N
    A = A / np.arange(N - 1, 0, -1)
    z = np.cumsum(np.hstack([zn, A]))[::-1]

    rho = np.diag(z / 2).astype(complex)

    rho[np.tril_indices(N, -1)] = X + 1j * Y
    rho += dagger(rho)

    return normalize(rho)


def Hermitian2v(H):
    N = H.shape[0]
    indices = np.triu_indices(N, 1)
    return np.hstack((H[indices].real, H[indices].imag, np.diag(H).real))


def v2Hermitian(V):
    N = int(round(np.sqrt(len(V))))

    X = V[:(N**2 - N) // 2]
    Y = V[(N**2 - N) // 2:(N**2 - N)]
    Z = V[(N**2 - N):]
    H = np.diag(Z / 2).astype(complex)

    H[np.triu_indices(N, 1)] = X + 1j * Y
    H += dagger(H)
    return H


def unitary2v(U):
    H = -1j * logm(U)
    return Hermitian2v(H)


def v2unitary(V):
    H = v2Hermitian(V)
    return expm(1j * H)


def rho2bloch(rho):
    """
    与 rho2v 类似，只不过是选择以 sigma_i \\otimes sigma_j ... 为基底来展开,
    只适用于N比特态,即 rho 的形状只能是 (2 ** N, 2 ** N)
    """
    bases = [s0, s1, s2, s3]
    N = int(np.log2(rho.shape[0]))
    return np.asarray([
        np.sum(rho * reduce(np.kron, s).T).real
        for s in product(bases, repeat=N)
    ][1:])


def bloch2rho(V):
    """
    rho2bloch 的逆函数
    """
    bases = [s0, s1, s2, s3]
    N = int(np.log2(len(V) + 1) / 2)
    rho = reduce(np.add,
                 (a * reduce(np.kron, s)
                  for (a, s) in zip(chain([1], V), product(bases, repeat=N))))
    return rho / rho.shape[0]


def traceDistance(rho: np.ndarray, sigma: np.ndarray) -> float:
    """Compute the trace distance between matrixes rho and sigma
    See Nielsen and Chuang, p. 403
    """
    A = rho - sigma
    return np.real(np.trace(sqrtm(dagger(A) @ A))) / 2.0


def fidelity(rho: np.ndarray, sigma: np.ndarray) -> float:
    '''
    The fidelity of two quantum states.
    '''
    rho = psi2rho(rho)
    sigma = psi2rho(sigma)
    rhosqrt = sqrtm(rho)
    return np.real(np.trace(sqrtm(rhosqrt @ sigma @ rhosqrt)))


def entropy(rho: np.ndarray) -> float:
    '''von Neumann / Shannon entropy function'''
    if rho.ndim == 1:  # pure states have no entropy
        return 0
    v = np.linalg.eigvalsh(rho).real
    return -np.sum(v[v > 0] * np.log2(v[v > 0]))


def concurrence(rho: np.ndarray) -> float:
    """Concurrence of a two-qubit density matrix.
    see http://qwiki.stanford.edu/wiki/Entanglement_of_Formation
    """
    yy = np.array([[ 0, 0, 0,-1],
                   [ 0, 0, 1, 0],
                   [ 0, 1, 0, 0],
                   [-1, 0, 0, 0]], dtype=complex) #yapf: disable
    m = rho @ yy @ rho.conj() @ yy
    eigs = [np.abs(e) for e in np.linalg.eig(m)[0]]
    e = [np.sqrt(x) for x in sorted(eigs, reverse=True)]
    return max(0, e[0] - e[1] - e[2] - e[3])


def eof(rho: np.ndarray) -> float:
    """Entanglement of formation of a two-qubit density matrix.
    see http://qwiki.stanford.edu/wiki/Entanglement_of_Formation
    """

    def h(x):
        if x <= 0 or x >= 1:
            return 0
        return -x * np.log2(x) - (1 - x) * np.log2(1 - x)

    C = concurrence(rho)
    arg = max(0, np.sqrt(1 - C**2))
    return h((1 + arg) / 2.0)


def randomUnitary(N: int) -> np.ndarray:
    """Random unitary matrix of dimension N."""
    H = np.random.randn(N, N) + 1j * np.random.randn(N, N)
    H = (H + dagger(H)) / 2
    U = expm(-1j * H)
    return U


def randomState(N: int) -> np.ndarray:
    """Generates a random pure state."""
    psi = np.random.randn(N) + 1j * np.random.randn(N)
    return normalize(psi)


def randomDensity(N: int, rank: int | None = None) -> np.ndarray:
    """Generates a random density matrix.  The distribution does not
    necessarily mean anything.  N is the dimension and rank is the
    number of non-zero eigenvalues."""

    if rank is None or rank > N:
        rank = N
    A = np.random.randn(N, rank) + 1j * np.random.randn(N, rank)
    rho = A @ A.T.conj()
    return rho / np.trace(rho)


##########################################################


def rabi(t, TR, Omega, A, offset):
    return oscillation(t, [(1, Omega)], A, offset) * decay(t, [TR])


def cpmg(t, T1, Tphi, Delta, A, offset, phi=0):
    return oscillation(t, [(np.exp(1j * phi), Delta)], A, offset) * decay(
        t, [2 * T1, Tphi])


def relaxation(t, T1, A, offset):
    return A * decay(t, [T1]) + offset

1	from functools import reduce	9✔
2	from itertools import chain, product	9✔
3
4	import numpy as np	9✔
5	from scipy.linalg import eigh, expm, logm, sqrtm	9✔
6
7	from waveforms.math.matricies import (BellPhiM, BellPhiP, BellPsiM, BellPsiP,	9✔
8	sigmaI, sigmaX, sigmaY, sigmaZ)
9	from waveforms.math.signal import decay, oscillation	9✔
10
11	# Paulis
12	s0, s1, s2, s3 = sigmaI(), sigmaX(), sigmaY(), sigmaZ()	9✔
13
14	# Bell states
15	phiplus, phiminus = BellPhiP(), BellPhiM()	9✔
16	psiplus, psiminus = BellPsiP(), BellPsiM()	9✔
17
18
19	def tensor(matrixes: np.ndarray) -> np.ndarray:	9✔
20	"""Compute the tensor product of a list (or array) of matrixes"""
UNCOV 21	return reduce(np.kron, matrixes)	×
22
23
24	def sigma(j: int, N: int = 1) -> np.ndarray:	9✔
25	"""
26	"""
UNCOV 27	s = [s0, s1, s2, s3]	×
UNCOV 28	dims = [4] * N	×
UNCOV 29	idx = np.unravel_index(j, dims)	×
UNCOV 30	return tensor(s[x] for x in idx)	×
31
32
33	def basis(x: int \| str, dim: int = 2) -> np.ndarray:	9✔
34	'''
35	Returns a single element of a state vector. Each component is
36	either a digit (0,1,2,...) representing the computational basis,
37	or one or two letters representing a special state (H,V,D,A,R,L,
38	X,Y,Z,Xm,Ym,Zm. B1..4 are the 4 bell states
39	'''
UNCOV 40	d = {	×
41	'H': [1, 0],
42	'V': [0, 1],
43	'D': [1, 1],
44	'A': [1, -1],
45	'R': [1, 1j],
46	'L': [1, -1j],
47	'X': [1, 1],
48	'Xp': [1, 1],
49	'Xm': [1, -1],
50	'Y': [1, 1j],
51	'Yp': [1, 1j],
52	'Ym': [1, -1j],
53	'Z': [1, 0],
54	'Zp': [1, 0],
55	'Zm': [0, 1],
56	'B1': phiminus,
57	'B2': phiplus,
58	'B3': psiminus,
59	'B4': psiplus
60	}
UNCOV 61	if isinstance(x, int) or x.isdigit():	×
UNCOV 62	rv = np.zeros((dim), dtype=complex)	×
UNCOV 63	rv[int(x)] = 1	×
64	else:
UNCOV 65	rv = np.array(d[x], dtype=complex)	×
66	rv = rv / np.linalg.norm(rv)	×
67	return rv	×
68
69
70	def commutator(A: np.ndarray, B: np.ndarray) -> np.ndarray:	9✔
71	"""Compute the commutator of two matrices"""
72	return A @ B - B @ A	×
73
74
75	def dagger(vector: np.ndarray) -> np.ndarray:	9✔
76	"""Compute the hermitian conjugate of a vector or matrix"""
77	return vector.conjugate().transpose()	9✔
78
79
80	def normalize(X: np.ndarray) -> np.ndarray:	9✔
81	"""Normalize the representation of X. This means:
82	For state vectors, X' * X = 1 and the first non-zero
83	element of X is positive real. For density matricies,
84	X is hermetian, unit trace, and non-negative"""
85	if X.ndim == 1:	9✔
UNCOV 86	nz = X[np.nonzero(X)][0]	×
UNCOV 87	X *= nz.conj()	×
UNCOV 88	amp = np.abs(dagger(X) @ X)	×
UNCOV 89	X /= np.sqrt(amp)	×
90	else:
91	(d, U) = eigh(X)	9✔
92	d = d.real	9✔
93	d[d < 0] = 0	9✔
94	X = U @ np.diag(d) @ dagger(U)	9✔
95	X = dagger(X) + X	9✔
96	X /= np.trace(X)	9✔
97	return X	9✔
98
99
100	def psi2rho(psi: np.ndarray) -> np.ndarray:	9✔
101	'''
102	Converts input state vector to density matrix represenation. If psi
103	is already a density matrix, leave it alone, so this function can be
104	used to coerce input arguments to the desired type.
105	'''
UNCOV 106	if psi.ndim == 1:	×
UNCOV 107	return normalize(np.outer(psi, psi.conj()))	×
108	else: # Already a density matrix
UNCOV 109	return normalize(psi)	×
110
111
112	def rho2v(rho):	9✔
113	"""密度矩阵转相干矢
114	"""
UNCOV 115	N = rho.shape[0]	×
UNCOV 116	indices = np.tril_indices(N, -1)	×
UNCOV 117	z = np.diag(rho).real	×
UNCOV 118	Z = np.cumsum(z[:-1]) - z[1:] * np.arange(1, N)	×
UNCOV 119	return np.hstack((rho[indices].real, rho[indices].imag, Z))	×
120
121
122	def v2rho(V):	9✔
123	"""相干矢转密度矩阵"""
124	N = int(np.sqrt(len(V) + 1))	×
UNCOV 125	assert N**2 - 1 == len(V)	×
126
UNCOV 127	X = V[:(N**2 - N) // 2]	×
UNCOV 128	Y = V[(N2 - N) // 2:(N2 - N)]	×
129	Z = V[(N**2 - N):]	×
130
UNCOV 131	A = np.hstack([Z[:0:-1] - Z[-2::-1], Z[0]])	×
132	zn = (1 - A.sum()) / N	×
133	A = A / np.arange(N - 1, 0, -1)	×
134	z = np.cumsum(np.hstack([zn, A]))[::-1]	×
135
136	rho = np.diag(z / 2).astype(complex)	×
137
138	rho[np.tril_indices(N, -1)] = X + 1j * Y	×
139	rho += dagger(rho)	×
140
141	return normalize(rho)	×
142
143
144	def Hermitian2v(H):	9✔
UNCOV 145	N = H.shape[0]	×
146	indices = np.triu_indices(N, 1)	×
UNCOV 147	return np.hstack((H[indices].real, H[indices].imag, np.diag(H).real))	×
148
149
150	def v2Hermitian(V):	9✔
151	N = int(round(np.sqrt(len(V))))	×
152
UNCOV 153	X = V[:(N**2 - N) // 2]	×
UNCOV 154	Y = V[(N2 - N) // 2:(N2 - N)]	×
UNCOV 155	Z = V[(N**2 - N):]	×
156	H = np.diag(Z / 2).astype(complex)	×
157
158	H[np.triu_indices(N, 1)] = X + 1j * Y	×
159	H += dagger(H)	×
160	return H	×
161
162
163	def unitary2v(U):	9✔
164	H = -1j * logm(U)	×
165	return Hermitian2v(H)	×
166
167
168	def v2unitary(V):	9✔
169	H = v2Hermitian(V)	×
170	return expm(1j * H)	×
171
172
173	def rho2bloch(rho):	9✔
174	"""
175	与 rho2v 类似，只不过是选择以 sigma_i \\otimes sigma_j ... 为基底来展开,
176	只适用于N比特态,即 rho 的形状只能是 (2 N, 2 N)
177	"""
UNCOV 178	bases = [s0, s1, s2, s3]	×
UNCOV 179	N = int(np.log2(rho.shape[0]))	×
UNCOV 180	return np.asarray([	×
181	np.sum(rho * reduce(np.kron, s).T).real
182	for s in product(bases, repeat=N)
183	][1:])
184
185
186	def bloch2rho(V):	9✔
187	"""
188	rho2bloch 的逆函数
189	"""
UNCOV 190	bases = [s0, s1, s2, s3]	×
UNCOV 191	N = int(np.log2(len(V) + 1) / 2)	×
UNCOV 192	rho = reduce(np.add,	×
193	(a * reduce(np.kron, s)
194	for (a, s) in zip(chain([1], V), product(bases, repeat=N))))
195	return rho / rho.shape[0]	×
196
197
198	def traceDistance(rho: np.ndarray, sigma: np.ndarray) -> float:	9✔
199	"""Compute the trace distance between matrixes rho and sigma
200	See Nielsen and Chuang, p. 403
201	"""
UNCOV 202	A = rho - sigma	×
UNCOV 203	return np.real(np.trace(sqrtm(dagger(A) @ A))) / 2.0	×
204
205
206	def fidelity(rho: np.ndarray, sigma: np.ndarray) -> float:	9✔
207	'''
208	The fidelity of two quantum states.
209	'''
UNCOV 210	rho = psi2rho(rho)	×
UNCOV 211	sigma = psi2rho(sigma)	×
UNCOV 212	rhosqrt = sqrtm(rho)	×
UNCOV 213	return np.real(np.trace(sqrtm(rhosqrt @ sigma @ rhosqrt)))	×
214
215
216	def entropy(rho: np.ndarray) -> float:	9✔
217	'''von Neumann / Shannon entropy function'''
218	if rho.ndim == 1: # pure states have no entropy	×
UNCOV 219	return 0	×
UNCOV 220	v = np.linalg.eigvalsh(rho).real	×
UNCOV 221	return -np.sum(v[v > 0] * np.log2(v[v > 0]))	×
222
223
224	def concurrence(rho: np.ndarray) -> float:	9✔
225	"""Concurrence of a two-qubit density matrix.
226	see http://qwiki.stanford.edu/wiki/Entanglement_of_Formation
227	"""
UNCOV 228	yy = np.array([[ 0, 0, 0,-1],	×
229	[ 0, 0, 1, 0],
230	[ 0, 1, 0, 0],
231	[-1, 0, 0, 0]], dtype=complex) #yapf: disable
UNCOV 232	m = rho @ yy @ rho.conj() @ yy	×
233	eigs = [np.abs(e) for e in np.linalg.eig(m)[0]]	×
UNCOV 234	e = [np.sqrt(x) for x in sorted(eigs, reverse=True)]	×
UNCOV 235	return max(0, e[0] - e[1] - e[2] - e[3])	×
236
237
238	def eof(rho: np.ndarray) -> float:	9✔
239	"""Entanglement of formation of a two-qubit density matrix.
240	see http://qwiki.stanford.edu/wiki/Entanglement_of_Formation
241	"""
242
UNCOV 243	def h(x):	×
UNCOV 244	if x <= 0 or x >= 1:	×
UNCOV 245	return 0	×
UNCOV 246	return -x * np.log2(x) - (1 - x) * np.log2(1 - x)	×
247
248	C = concurrence(rho)	×
249	arg = max(0, np.sqrt(1 - C**2))	×
250	return h((1 + arg) / 2.0)	×
251
252
253	def randomUnitary(N: int) -> np.ndarray:	9✔
254	"""Random unitary matrix of dimension N."""
255	H = np.random.randn(N, N) + 1j * np.random.randn(N, N)	9✔
256	H = (H + dagger(H)) / 2	9✔
257	U = expm(-1j * H)	9✔
258	return U	9✔
259
260
261	def randomState(N: int) -> np.ndarray:	9✔
262	"""Generates a random pure state."""
UNCOV 263	psi = np.random.randn(N) + 1j * np.random.randn(N)	×
UNCOV 264	return normalize(psi)	×
265
266
267	def randomDensity(N: int, rank: int \| None = None) -> np.ndarray:	9✔
268	"""Generates a random density matrix. The distribution does not
269	necessarily mean anything. N is the dimension and rank is the
270	number of non-zero eigenvalues."""
271
272	if rank is None or rank > N:	9✔
273	rank = N	9✔
274	A = np.random.randn(N, rank) + 1j * np.random.randn(N, rank)	9✔
275	rho = A @ A.T.conj()	9✔
276	return rho / np.trace(rho)	9✔
277
278
279	##########################################################
280
281
282	def rabi(t, TR, Omega, A, offset):	9✔
UNCOV 283	return oscillation(t, [(1, Omega)], A, offset) * decay(t, [TR])	×
284
285
286	def cpmg(t, T1, Tphi, Delta, A, offset, phi=0):	9✔
UNCOV 287	return oscillation(t, [(np.exp(1j * phi), Delta)], A, offset) * decay(	×
288	t, [2 * T1, Tphi])
289
290
291	def relaxation(t, T1, A, offset):	9✔
292	return A * decay(t, [T1]) + offset	×

feihoo87 / waveforms / 6893865468

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