6534953321

Committed 16 Oct 2023 02:19PM UTC coverage: 35.674% (-22.7%) from 58.421%

Build # 6534953321

Build Type

push

github

Committed by

feihoo87

Commit Message

fix Coveralls

Run Details

5913 of 16575 relevant lines covered (35.67%)

3.21 hits per line

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

29.79

/waveforms/qlisp/simulator/__init__.py

import itertools
from functools import partial, reduce

import numpy as np

from .mat import U, fSim, make_immutable, rfUnitary

try:
    import tensornetwork as tn
except ImportError:
    pass

__matrix_of_gates = {}


def regesterGateMatrix(gate, mat, N=None):
    if isinstance(mat, np.ndarray):
        mat = make_immutable(mat)
    if N is None:
        N = int(np.log2(mat.shape[0]))
    __matrix_of_gates[gate] = (mat, N)


def gate2mat(gate):
    if isinstance(gate, str) and gate in __matrix_of_gates:
        return __matrix_of_gates[gate][0]
    elif isinstance(gate, tuple) and gate[0] in __matrix_of_gates:
        if callable(__matrix_of_gates[gate[0]][0]):
            return __matrix_of_gates[gate[0]][0](*gate[1:])
        else:
            raise ValueError(
                f"Could not call {gate[0]}(*{gate[1:]}), `{gate[0]}` is not callable."
            )
    elif isinstance(gate, tuple) and gate[0] == 'C':
        U = gate2mat(gate[1])
        ret = np.eye(2 * U.shape[0], dtype=np.complex)
        ret[U.shape[0]:, U.shape[0]:] = U
        return ret
    else:
        raise ValueError(f'Unexcept gate {gate}')


def splite_at(l, bits):
    """将 l 的二进制位于 bits 所列的位置上断开插上0
    
    如 splite_at(int('1111',2), [0,2,4,6]) == int('10101010', 2)
    bits 必须从小到大排列
    """
    r = l
    for n in bits:
        mask = (1 << n) - 1
        low = r & mask
        high = r - low
        r = (high << 1) + low
    return r


def place_at(l, bits):
    """将 l 的二进制位置于 bits 所列的位置上
    
    如 place_at(int('10111',2), [0,2,4,5,6]) == int('01010101', 2)
    """
    r = 0
    for index, n in enumerate(bits):
        b = (l >> index) & 1
        r += b << n
    return r


def reduceSubspace(targets, N, inputMat, func, args):
    innerDim = 2**len(targets)
    outerDim = 2**(N - len(targets))

    targets = tuple(reversed([N - i - 1 for i in targets]))

    def index(targets, i, j):
        return splite_at(j, sorted(targets)) | place_at(i, targets)

    if len(inputMat.shape) == 1:
        for k in range(outerDim):
            innerIndex = [index(targets, i, k) for i in range(innerDim)]
            inputMat[innerIndex] = func(inputMat[innerIndex], *args)
    else:
        for k, l in itertools.product(range(outerDim), repeat=2):
            innerIndex = np.asarray(
                [[index(targets, i, k),
                  index(targets, j, l)]
                 for i, j in itertools.product(range(innerDim), repeat=2)]).T
            sub = inputMat[innerIndex[0], innerIndex[1]].reshape(
                (innerDim, innerDim))
            inputMat[innerIndex[0], innerIndex[1]] = func(sub, *args).flatten()
    return inputMat


def applySeq(seq, psi0=None):
    if psi0 is None:
        psi = np.array([1, 0], dtype=complex)
        N = 1
    else:
        psi = psi0
        N = int(np.round(np.log2(psi.shape[0])))

    if len(psi.shape) == 1:
        I = np.array([1, 0])
    else:
        I = np.eye(2)

    def func(psi, U):
        return U @ psi

    for gate, qubits in seq:
        if (isinstance(gate, tuple) and gate[0] in ['Measure', 'Delay']) or (
                isinstance(gate, str) and gate in ['Barrier']):
            continue
        if isinstance(qubits, tuple):
            M = max(qubits)
        else:
            M = qubits
            qubits = (qubits, )
        if M >= N:
            psi = reduce(np.kron, itertools.repeat(I, times=M - N + 1), psi)
            N = M + 1

        reduceSubspace(qubits, N, psi, func, (gate2mat(gate), ))

    return psi


def seq2mat(seq):
    return applySeq(seq, np.eye(2, dtype=complex))


def apply_gate(qubit_edges, gate, operating_qubits):
    if 'tn' not in globals():
        raise ImportError('Please install tensornetwork first.')

    op = tn.Node(gate)
    for i, bit in enumerate(operating_qubits):
        tn.connect(qubit_edges[bit], op[i])
        qubit_edges[bit] = op[i + len(operating_qubits)]


def circuit_network(circ):
    if 'tn' not in globals():
        raise ImportError('Please install tensornetwork first.')

    N = 0

    all_nodes = []

    with tn.NodeCollection(all_nodes):
        left_edges = []
        right_edges = []

        for gate, qubits in circ:
            if (isinstance(gate, tuple)
                    and gate[0] in ['Measure', 'Delay']) or (isinstance(
                        gate, str) and gate in ['Barrier']):
                continue
            if isinstance(qubits, tuple):
                M = max(qubits)
            else:
                M = qubits
                qubits = (qubits, )
            if M >= N:
                new_nodes = [tn.Node(np.eye(2)) for _ in range(M - N + 1)]
                left_edges.extend([node[0] for node in new_nodes])
                right_edges.extend([node[1] for node in new_nodes])
                N = M + 1

            gate_mat = gate2mat(gate)
            gate_tenser = gate_mat.reshape((2, 2) * len(qubits))
            apply_gate(left_edges, gate_tenser, qubits)

    return all_nodes, left_edges, right_edges, N


def apply_circuit(circ, qubits=None, init_state=None):
    if 'tn' not in globals():
        raise ImportError('Please install tensornetwork first.')

    all_nodes, left_edges, right_edges, N = circuit_network(circ)
    if qubits is not None:
        assert len(qubits) == N
    else:
        qubits = list(range(N))
    if init_state is None:
        init_state_nodes = [
            tn.Node(np.array([1, 0], dtype=complex)) for i in range(N)
        ]
        init_state_edges = [init_state_nodes[i][0] for i in range(N)]
    else:
        init_state_nodes, init_state_edges = init_state

    all_nodes.extend(init_state_nodes)
    qubit_edges = [left_edges[q] for q in qubits]
    for i, q in enumerate(qubits):
        tn.connect(right_edges[i], init_state_edges[q])
    result = tn.contractors.optimal(all_nodes, output_edge_order=qubit_edges)
    return result.tensor.reshape(-1)


def circuit2mat(circ):
    if 'tn' not in globals():
        raise ImportError('Please install tensornetwork first.')

    all_nodes, left_edges, right_edges, N = circuit_network(circ)
    result = tn.contractors.optimal(all_nodes,
                                    output_edge_order=left_edges + right_edges)
    return result.tensor.reshape((2**N, 2**N))


regesterGateMatrix('U', U, 1)
regesterGateMatrix('P', lambda p: U(theta=0, phi=0, lambda_=p), 1)
regesterGateMatrix('rfUnitary', rfUnitary, 1)
regesterGateMatrix('Rx', partial(rfUnitary, phi=0), 1)
regesterGateMatrix('Ry', partial(rfUnitary, phi=np.pi / 2), 1)
regesterGateMatrix('Rz', lambda p: U(theta=0, phi=0, lambda_=p), 1)
regesterGateMatrix('fSim', fSim, 2)

# one qubit
regesterGateMatrix('I', np.array([[1, 0], [0, 1]]))
regesterGateMatrix('X', np.array([[0, -1j], [-1j, 0]]))
regesterGateMatrix('Y', np.array([[0, -1], [1, 0]]))
regesterGateMatrix('X/2', np.array([[1, -1j], [-1j, 1]]) / np.sqrt(2))
regesterGateMatrix('Y/2', np.array([[1, -1], [1, 1]]) / np.sqrt(2))
regesterGateMatrix('-X/2', np.array([[1, 1j], [1j, 1]]) / np.sqrt(2))
regesterGateMatrix('-Y/2', np.array([[1, 1], [-1, 1]]) / np.sqrt(2))
regesterGateMatrix('Z', np.array([[1, 0], [0, -1]]))
regesterGateMatrix('S', np.array([[1, 0], [0, 1j]]))
regesterGateMatrix('-S', np.array([[1, 0], [0, -1j]]))
regesterGateMatrix('H', np.array([[1, 1], [1, -1]]) / np.sqrt(2))

# non-clifford
regesterGateMatrix('T',
                   np.array([[1, 0], [0, 1 / np.sqrt(2) + 1j / np.sqrt(2)]]))
regesterGateMatrix('-T',
                   np.array([[1, 0], [0, 1 / np.sqrt(2) - 1j / np.sqrt(2)]]))
regesterGateMatrix('W/2', rfUnitary(np.pi / 2, np.pi / 4))
regesterGateMatrix('-W/2', rfUnitary(-np.pi / 2, np.pi / 4))
regesterGateMatrix('V/2', rfUnitary(np.pi / 2, 3 * np.pi / 4))
regesterGateMatrix('-V/2', rfUnitary(-np.pi / 2, 3 * np.pi / 4))

# two qubits
regesterGateMatrix(
    'CZ', np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, -1]]))
regesterGateMatrix(
    'Cnot', np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]]))
regesterGateMatrix(
    'iSWAP',
    np.array([[1, 0, 0, 0], [0, 0, 1j, 0], [0, 1j, 0, 0], [0, 0, 0, 1]]))
regesterGateMatrix(
    'SWAP', np.array([[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]]))
regesterGateMatrix(
    'CR',
    np.array([[1, 1j, 0, 0], [1j, 1, 0, 0], [0, 0, 1, -1j], [0, 0, -1j, 1]]) /
    np.sqrt(2))

# non-clifford
regesterGateMatrix(
    'SQiSWAP',
    np.array([[1, 0, 0, 0], [0, 1 / np.sqrt(2), 1j / np.sqrt(2), 0],
              [0, 1j / np.sqrt(2), 1 / np.sqrt(2), 0], [0, 0, 0, 1]]))

if __name__ == '__main__':
    # Porter-Thomas distribution

    def randomSeq(depth, N):
        seq = []
        for i in range(depth):
            for j in range(N):
                seq.append((np.random.choice(['X/2', 'Y/2', 'W/2']), j))
            for j in range(i % 2, N, 2):
                seq.append(('SQiSWAP', (j, (j + 1) % N)))
        return seq

    p = []
    # run 1000 random circuit on 6 qubits
    for i in range(1000):
        print('    ', i, end='')
        seq = randomSeq(50, 6)
        psi = applySeq(seq)
        p.extend(list(np.abs(psi)**2))
        print('    ', i)
    p = np.asarray(p)

    # plot distribution of probabilities
    N = 2**6
    y, x = np.histogram(N * p, bins=50, density=True)

    import matplotlib.pyplot as plt

    plt.semilogy((x[:-1] + x[1:]) / 2, y)
    plt.show()

1	import itertools	9✔
2	from functools import partial, reduce	9✔
3
4	import numpy as np	9✔
5
6	from .mat import U, fSim, make_immutable, rfUnitary	9✔
7
8	try:	9✔
9	import tensornetwork as tn	9✔
10	except ImportError:	9✔
11	pass	9✔
12
13	__matrix_of_gates = {}	9✔
14
15
16	def regesterGateMatrix(gate, mat, N=None):	9✔
17	if isinstance(mat, np.ndarray):	9✔
18	mat = make_immutable(mat)	9✔
19	if N is None:	9✔
20	N = int(np.log2(mat.shape[0]))	9✔
21	__matrix_of_gates[gate] = (mat, N)	9✔
22
23
24	def gate2mat(gate):	9✔
25	if isinstance(gate, str) and gate in __matrix_of_gates:	×
26	return __matrix_of_gates[gate][0]	×
27	elif isinstance(gate, tuple) and gate[0] in __matrix_of_gates:	×
28	if callable(__matrix_of_gates[gate[0]][0]):	×
29	return __matrix_of_gates[gate[0]][0](*gate[1:])	×
30	else:
31	raise ValueError(	×
32	f"Could not call {gate[0]}(*{gate[1:]}), `{gate[0]}` is not callable."
33	)
34	elif isinstance(gate, tuple) and gate[0] == 'C':	×
35	U = gate2mat(gate[1])	×
36	ret = np.eye(2 * U.shape[0], dtype=np.complex)	×
37	ret[U.shape[0]:, U.shape[0]:] = U	×
38	return ret	×
39	else:
40	raise ValueError(f'Unexcept gate {gate}')	×
41
42
43	def splite_at(l, bits):	9✔
44	"""将 l 的二进制位于 bits 所列的位置上断开插上0
45
46	如 splite_at(int('1111',2), [0,2,4,6]) == int('10101010', 2)
47	bits 必须从小到大排列
48	"""
49	r = l	×
50	for n in bits:	×
51	mask = (1 << n) - 1	×
52	low = r & mask	×
53	high = r - low	×
54	r = (high << 1) + low	×
55	return r	×
56
57
58	def place_at(l, bits):	9✔
59	"""将 l 的二进制位置于 bits 所列的位置上
60
61	如 place_at(int('10111',2), [0,2,4,5,6]) == int('01010101', 2)
62	"""
63	r = 0	×
64	for index, n in enumerate(bits):	×
65	b = (l >> index) & 1	×
66	r += b << n	×
67	return r	×
68
69
70	def reduceSubspace(targets, N, inputMat, func, args):	9✔
71	innerDim = 2**len(targets)	×
72	outerDim = 2**(N - len(targets))	×
73
74	targets = tuple(reversed([N - i - 1 for i in targets]))	×
75
76	def index(targets, i, j):	×
77	return splite_at(j, sorted(targets)) \| place_at(i, targets)	×
78
79	if len(inputMat.shape) == 1:	×
80	for k in range(outerDim):	×
81	innerIndex = [index(targets, i, k) for i in range(innerDim)]	×
82	inputMat[innerIndex] = func(inputMat[innerIndex], *args)	×
83	else:
84	for k, l in itertools.product(range(outerDim), repeat=2):	×
85	innerIndex = np.asarray(	×
86	[[index(targets, i, k),
87	index(targets, j, l)]
88	for i, j in itertools.product(range(innerDim), repeat=2)]).T
89	sub = inputMat[innerIndex[0], innerIndex[1]].reshape(	×
90	(innerDim, innerDim))
91	inputMat[innerIndex[0], innerIndex[1]] = func(sub, *args).flatten()	×
92	return inputMat	×
93
94
95	def applySeq(seq, psi0=None):	9✔
96	if psi0 is None:	×
97	psi = np.array([1, 0], dtype=complex)	×
98	N = 1	×
99	else:
100	psi = psi0	×
101	N = int(np.round(np.log2(psi.shape[0])))	×
102
103	if len(psi.shape) == 1:	×
104	I = np.array([1, 0])	×
105	else:
106	I = np.eye(2)	×
107
108	def func(psi, U):	×
109	return U @ psi	×
110
111	for gate, qubits in seq:	×
112	if (isinstance(gate, tuple) and gate[0] in ['Measure', 'Delay']) or (	×
113	isinstance(gate, str) and gate in ['Barrier']):
114	continue	×
115	if isinstance(qubits, tuple):	×
116	M = max(qubits)	×
117	else:
118	M = qubits	×
119	qubits = (qubits, )	×
120	if M >= N:	×
121	psi = reduce(np.kron, itertools.repeat(I, times=M - N + 1), psi)	×
122	N = M + 1	×
123
124	reduceSubspace(qubits, N, psi, func, (gate2mat(gate), ))	×
125
126	return psi	×
127
128
129	def seq2mat(seq):	9✔
130	return applySeq(seq, np.eye(2, dtype=complex))	×
131
132
133	def apply_gate(qubit_edges, gate, operating_qubits):	9✔
134	if 'tn' not in globals():	×
135	raise ImportError('Please install tensornetwork first.')	×
136
137	op = tn.Node(gate)	×
138	for i, bit in enumerate(operating_qubits):	×
139	tn.connect(qubit_edges[bit], op[i])	×
140	qubit_edges[bit] = op[i + len(operating_qubits)]	×
141
142
143	def circuit_network(circ):	9✔
144	if 'tn' not in globals():	×
145	raise ImportError('Please install tensornetwork first.')	×
146
147	N = 0	×
148
149	all_nodes = []	×
150
151	with tn.NodeCollection(all_nodes):	×
152	left_edges = []	×
153	right_edges = []	×
154
155	for gate, qubits in circ:	×
156	if (isinstance(gate, tuple)	×
157	and gate[0] in ['Measure', 'Delay']) or (isinstance(
158	gate, str) and gate in ['Barrier']):
159	continue	×
160	if isinstance(qubits, tuple):	×
161	M = max(qubits)	×
162	else:
163	M = qubits	×
164	qubits = (qubits, )	×
165	if M >= N:	×
166	new_nodes = [tn.Node(np.eye(2)) for _ in range(M - N + 1)]	×
167	left_edges.extend([node[0] for node in new_nodes])	×
168	right_edges.extend([node[1] for node in new_nodes])	×
169	N = M + 1	×
170
171	gate_mat = gate2mat(gate)	×
172	gate_tenser = gate_mat.reshape((2, 2) * len(qubits))	×
173	apply_gate(left_edges, gate_tenser, qubits)	×
174
175	return all_nodes, left_edges, right_edges, N	×
176
177
178	def apply_circuit(circ, qubits=None, init_state=None):	9✔
179	if 'tn' not in globals():	×
180	raise ImportError('Please install tensornetwork first.')	×
181
182	all_nodes, left_edges, right_edges, N = circuit_network(circ)	×
183	if qubits is not None:	×
184	assert len(qubits) == N	×
185	else:
186	qubits = list(range(N))	×
187	if init_state is None:	×
188	init_state_nodes = [	×
189	tn.Node(np.array([1, 0], dtype=complex)) for i in range(N)
190	]
191	init_state_edges = [init_state_nodes[i][0] for i in range(N)]	×
192	else:
193	init_state_nodes, init_state_edges = init_state	×
194
195	all_nodes.extend(init_state_nodes)	×
196	qubit_edges = [left_edges[q] for q in qubits]	×
197	for i, q in enumerate(qubits):	×
198	tn.connect(right_edges[i], init_state_edges[q])	×
199	result = tn.contractors.optimal(all_nodes, output_edge_order=qubit_edges)	×
200	return result.tensor.reshape(-1)	×
201
202
203	def circuit2mat(circ):	9✔
204	if 'tn' not in globals():	×
205	raise ImportError('Please install tensornetwork first.')	×
206
207	all_nodes, left_edges, right_edges, N = circuit_network(circ)	×
208	result = tn.contractors.optimal(all_nodes,	×
209	output_edge_order=left_edges + right_edges)
210	return result.tensor.reshape((2N, 2N))	×
211
212
213	regesterGateMatrix('U', U, 1)	9✔
214	regesterGateMatrix('P', lambda p: U(theta=0, phi=0, lambda_=p), 1)	9✔
215	regesterGateMatrix('rfUnitary', rfUnitary, 1)	9✔
216	regesterGateMatrix('Rx', partial(rfUnitary, phi=0), 1)	9✔
217	regesterGateMatrix('Ry', partial(rfUnitary, phi=np.pi / 2), 1)	9✔
218	regesterGateMatrix('Rz', lambda p: U(theta=0, phi=0, lambda_=p), 1)	9✔
219	regesterGateMatrix('fSim', fSim, 2)	9✔
220
221	# one qubit
222	regesterGateMatrix('I', np.array([[1, 0], [0, 1]]))	9✔
223	regesterGateMatrix('X', np.array([[0, -1j], [-1j, 0]]))	9✔
224	regesterGateMatrix('Y', np.array([[0, -1], [1, 0]]))	9✔
225	regesterGateMatrix('X/2', np.array([[1, -1j], [-1j, 1]]) / np.sqrt(2))	9✔
226	regesterGateMatrix('Y/2', np.array([[1, -1], [1, 1]]) / np.sqrt(2))	9✔
227	regesterGateMatrix('-X/2', np.array([[1, 1j], [1j, 1]]) / np.sqrt(2))	9✔
228	regesterGateMatrix('-Y/2', np.array([[1, 1], [-1, 1]]) / np.sqrt(2))	9✔
229	regesterGateMatrix('Z', np.array([[1, 0], [0, -1]]))	9✔
230	regesterGateMatrix('S', np.array([[1, 0], [0, 1j]]))	9✔
231	regesterGateMatrix('-S', np.array([[1, 0], [0, -1j]]))	9✔
232	regesterGateMatrix('H', np.array([[1, 1], [1, -1]]) / np.sqrt(2))	9✔
233
234	# non-clifford
235	regesterGateMatrix('T',	9✔
236	np.array([[1, 0], [0, 1 / np.sqrt(2) + 1j / np.sqrt(2)]]))
237	regesterGateMatrix('-T',	9✔
238	np.array([[1, 0], [0, 1 / np.sqrt(2) - 1j / np.sqrt(2)]]))
239	regesterGateMatrix('W/2', rfUnitary(np.pi / 2, np.pi / 4))	9✔
240	regesterGateMatrix('-W/2', rfUnitary(-np.pi / 2, np.pi / 4))	9✔
241	regesterGateMatrix('V/2', rfUnitary(np.pi / 2, 3 * np.pi / 4))	9✔
242	regesterGateMatrix('-V/2', rfUnitary(-np.pi / 2, 3 * np.pi / 4))	9✔
243
244	# two qubits
245	regesterGateMatrix(	9✔
246	'CZ', np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, -1]]))
247	regesterGateMatrix(	9✔
248	'Cnot', np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]]))
249	regesterGateMatrix(	9✔
250	'iSWAP',
251	np.array([[1, 0, 0, 0], [0, 0, 1j, 0], [0, 1j, 0, 0], [0, 0, 0, 1]]))
252	regesterGateMatrix(	9✔
253	'SWAP', np.array([[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]]))
254	regesterGateMatrix(	9✔
255	'CR',
256	np.array([[1, 1j, 0, 0], [1j, 1, 0, 0], [0, 0, 1, -1j], [0, 0, -1j, 1]]) /
257	np.sqrt(2))
258
259	# non-clifford
260	regesterGateMatrix(	9✔
261	'SQiSWAP',
262	np.array([[1, 0, 0, 0], [0, 1 / np.sqrt(2), 1j / np.sqrt(2), 0],
263	[0, 1j / np.sqrt(2), 1 / np.sqrt(2), 0], [0, 0, 0, 1]]))
264
265	if __name__ == '__main__':	9✔
266	# Porter-Thomas distribution
267
268	def randomSeq(depth, N):	×
269	seq = []	×
270	for i in range(depth):	×
271	for j in range(N):	×
272	seq.append((np.random.choice(['X/2', 'Y/2', 'W/2']), j))	×
273	for j in range(i % 2, N, 2):	×
274	seq.append(('SQiSWAP', (j, (j + 1) % N)))	×
275	return seq	×
276
277	p = []	×
278	# run 1000 random circuit on 6 qubits
279	for i in range(1000):	×
280	print(' ', i, end='')	×
281	seq = randomSeq(50, 6)	×
282	psi = applySeq(seq)	×
283	p.extend(list(np.abs(psi)**2))	×
284	print(' ', i)	×
285	p = np.asarray(p)	×
286
287	# plot distribution of probabilities
288	N = 2**6	×
289	y, x = np.histogram(N * p, bins=50, density=True)	×
290
291	import matplotlib.pyplot as plt	×
292
293	plt.semilogy((x[:-1] + x[1:]) / 2, y)	×
294	plt.show()	×

feihoo87 / waveforms / 6534953321

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