-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathvqcd_main_funcs.py
251 lines (198 loc) · 8.08 KB
/
vqcd_main_funcs.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
from extensions import *
from vqcd_secondary_funcs import *
# Cost function evaluator function
# Based Fig. 1a & 1b of: https://www.nature.com/articles/s41534-019-0167-6
def cost_func(qdim, theta, kraus_chan, ansatz_no, shots, purity_before_diag):
""" Calculate the value of the cost function for particular choice of the
input:
------
qdim = dimention of quantum channel.
theta = random initial angles.
channel = kraus operator for the quantum channel.
ansatz_no = ansatz number.
shots = number of shors for simulator.
purity_before_diag = purity of channel before diagonalization.
return:
-------
cost function value
"""
# number of qubits in the full circuit
n = 2 * qdim
# n is also the number of parameters in the ansatz
theta = theta.reshape(-1, 2 * n)
# start the ciruit
qc_cost = QuantumCircuit(2 * n, 2 * qdim)
# add the selected ansatz twice
qc_cost = qc_cost.compose(jamilchoi(qdim, kraus_chan), list(range(n)))
qc_cost = qc_cost.compose(jamilchoi(qdim, kraus_chan), list(range(n, 2 * n)))
if qdim == 1:
qc_cost = qc_cost.compose(ansatz_2q(n, theta, ansatz_no), list(range(n)))
qc_cost = qc_cost.compose(ansatz_2q(n, theta, ansatz_no), list(range(n, 2 * n)))
elif qdim == 2:
qc_cost = qc_cost.compose(ansatz_4q(n, theta, ansatz_no), list(range(n)))
qc_cost = qc_cost.compose(ansatz_4q(n, theta, ansatz_no), list(range(n, 2 * n)))
# add cx operations
for i in range(n):
qc_cost.cx([i + n], [i])
# evolve in input matrix
qc_cost.measure(list(range(2 * qdim)), list(range(2 * qdim)))
device = Aer.get_backend('qasm_simulator')
result = execute(qc_cost, backend=device, shots=shots).result()
counts = result.get_counts(qc_cost)
purity_after_diag = counts['0' * 2 * qdim] / shots
cost = purity_before_diag - purity_after_diag
return cost.real
# Estimated eigenvalues
# Based on Fig.1c of: https://www.nature.com/articles/s41534-019-0167-6
def eig_info(qdim, theta, kraus_chan, ansatz_no, error, device_type, noise_mdl, noise_amp):
"""
Return the measurement basis based on the relative error on calculating eigenvalues.
input:
------
qdim = dimension of quantum channel.
theta = optimized angles obtained corresponding to optimal cost.
kraus_chan = kraus operators for the quantum channel.
ansatz_no = the no. of ansatz under use.
error = threshold for relative error.
device_type = for simulated -->> 'sim,
for real -->> 'real.
noise_mdl =
***for 'sim'***, specify noise model,
amplitude damping -->> 'amp_damp',
depolarizing -->> 'depol',
random X -->> 'rand_x';
***for 'real'***, specify ibmq device,
e.g. 'ibmq_belem', 'ibmq_lima' etc.
noise_amp == amplitude of noise, 0 <= noise_amp <= 1.
#######################
output:
-------
a list containing inferred eigenvalues
"""
n = 2 * qdim
qc_meas = QuantumCircuit(n)
qc_meas = qc_meas.compose(jamilchoi(qdim, kraus_chan), list(range(n)))
if qdim == 1:
qc_meas = qc_meas.compose(ansatz_2q(n, theta, ansatz_no), list(range(n)))
elif qdim == 2:
qc_meas = qc_meas.compose(ansatz_4q(n, theta, ansatz_no), list(range(n)))
for q in range(n // 2):
qc_meas.swap(q, n - q - 1)
qc_meas.measure_all()
shots = 20000
if device_type == 'sim':
device = Aer.get_backend('qasm_simulator')
if noise_mdl == 'amp_damp':
result = execute(qc_meas, backend=device, noise_model=create_noise_model_amp_damping(noise_amp),
shots=shots).result()
elif noise_mdl == 'depol':
result = execute(qc_meas, backend=device, noise_model=create_noise_model_depolarazing(noise_amp),
shots=shots).result()
elif noise_mdl == 'rand_x':
result = execute(qc_meas, backend=device, noise_model=create_noise_model_random_x(noise_amp),
shots=shots).result()
elif noise_mdl == 'simulator':
result = execute(qc_meas, backend=device, shots=shots).result()
elif device_type == 'real':
shots = 2000
IBMQ.load_account()
provider = IBMQ.get_provider(hub='ibm-q')
device = provider.get_backend(noise_mdl)
sim_real = AerSimulator.from_backend(device)
result = sim_real.run(qc_meas).result()
counts = result.get_counts(qc_meas)
tmp_counts = {}
for num in range(0, 2 ** n):
tmp_counts[format(num, '0{}b'.format(n))] = 0
for ks in counts.keys():
tmp_counts[ks] = counts[ks]
sim_eigval = {}
for num in range(0, 2 ** n):
stat_count = tmp_counts[format(num, '0{}b'.format(n))]
rel_error = np.sqrt(shots) / stat_count
if rel_error <= error:
est_eig = stat_count / shots
sim_eigval[format(num, '0{}b'.format(n))] = est_eig
else:
sim_eigval[f'n{num}'] = 0
return sim_eigval
# Truncated Fidelity estimator function
# Described in Fig.1 of: https://arxiv.org/abs/1906.09253
def trun_output(n, any_state, channel, theta, ansatz_no, error, device_type, noise_mdl, noise_amp):
"""
Return the measurement basis based on the relative error on calculating eigenvalues.
input:
------
qdim = dimention of quantum channel.
theta = optimized angles obtained corresponding to optimal cost.
kraus_chan = kraus operators for the quantum channel.
ansatz_no = the no. of ansatz under use.
error = threshold for relative error.
device_type = for simulated -->> 'sim,
for real -->> 'real.
noise_mdl = for 'sim', specify noise model, amplitude damping -->> 'amp_damp', depolarizing -->> 'depol',
random X -->> 'rand_x'.
noise_amp = amplitude of noise, 0 <= noise_amp <= 1.
#######################
output:
-------
truncated fidelity counds i.e. TFB, TGFB values
"""
eig_dict = eig_info(n // 2, theta, channel, ansatz_no, error, device_type, noise_mdl, noise_amp)
basis = []
eigval = []
for key in eig_dict.keys():
basis.append(key)
eigval.append(eig_dict[key])
eigenvector = []
for i in basis:
if i == 'n':
eigenvector.append(np.zeros(2 ** n))
else:
qc_reg = QuantumRegister(n)
qc = QuantumCircuit(qc_reg)
for q in range(len(i)):
if i[q] == '0':
qc.id([q])
else:
qc.x([q])
if n // 2 == 1:
qc = qc.compose(ansatz_2q(n, theta, ansatz_no).inverse())
if n // 2 == 2:
qc = qc.compose(ansatz_4q(n, theta, ansatz_no).inverse())
simulator = Aer.get_backend('statevector_simulator')
result = execute(qc, simulator).result()
eigen_vec = result.get_statevector(qc)
eigenvector.append(eigen_vec)
T = np.zeros((2 ** n, 2 ** n))
Tlist = []
for r_i in range(len(eigval)):
for r_j in range(len(eigval)):
outprod = outerproduct(eigenvector[r_i], eigenvector[r_j])
expval = np.sum(np.diag(np.matmul(any_state, outprod)))
sqrt_eig = np.sqrt(eigval[r_i] * eigval[r_j])
Tij = np.multiply(sqrt_eig, expval)
Tlist.append(Tij)
T = np.add(T, np.multiply(Tij, outprod))
sigma = np.asarray(Tlist).reshape(-1, 2 ** n)
T = np.asmatrix(T)
eigT = la.eig(T)[0].real
flag = 0
for e in eigT:
if e < 0:
flag += 1
if flag != 0:
eigT, _ = maximum_likelihood(T)
# TFB calculation
TFB = np.sum(np.sqrt(np.abs(eigT)))
# TGFB calculation
if error == 1:
TGFB = TFB
else:
TGFB = TFB + np.sqrt((1 - sum(eigval)) * (1 - sum(np.diag(sigma))))
# counts rank corresponding to threshold of relative error
m = 0
for e in eigval:
if e != 0:
m += 1
return np.abs(TFB), np.abs(TGFB), m