-
Notifications
You must be signed in to change notification settings - Fork 9
/
dynamics.py
256 lines (190 loc) · 9.51 KB
/
dynamics.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
252
253
254
255
256
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
def make_heisenberg_circuits(n_spins, M, omegadt):
circuits = []
circuit = QuantumCircuit(n_spins)
# 第0ビットを 1/√2 (|0> + |1>) にする
circuit.h(0)
# Δtでの時間発展をM回繰り返すループ
for istep in range(M):
# ハミルトニアンのn-1個の項への分解に関するループ
for jspin in range(n_spins - 1):
# ZZ
circuit.cx(jspin, jspin + 1)
circuit.rz(-omegadt, jspin + 1)
circuit.cx(jspin, jspin + 1)
# XX YY
circuit.h([jspin, jspin + 1])
circuit.s([jspin, jspin + 1])
circuit.cx(jspin, jspin + 1)
circuit.rx(-omegadt, jspin)
circuit.rz(-omegadt, jspin + 1)
circuit.cx(jspin, jspin + 1)
circuit.sdg([jspin, jspin + 1])
circuit.h([jspin, jspin + 1])
# この時点での回路のコピーをリストに保存
# measure_all(inplace=False) はここまでの回路のコピーに測定を足したものを返す
circuits.append(circuit.measure_all(inplace=False))
return circuits
def bit_expectations_sv(time_points, statevectors):
"""Compute the bit expectation values at each time point from statevectors.
Args:
time_points (np.ndarray(shape=(T,), dtype=float)): Time points.
statevectors (np.ndarray(shape=(D, T), dtype=np.complex128)): State vector as a function of time.
Returns:
np.ndarray(shape=(T, n), dtype=float): Time points tiled for each bit.
np.ndarray(shape=(T, n), dtype=float): Bit expectation values.
"""
num_bits = np.round(np.log2(statevectors.shape[0])).astype(int)
if num_bits > 8:
raise NotImplementedError('Function not compatible with number of qubits > 8')
# Probability of seeing each bitstring at each time point
probs = np.square(np.abs(statevectors)) # shape (D, T)
# Unpack each index into a binary
indices = np.expand_dims(np.arange(2 ** num_bits, dtype=np.uint8), axis=1) # shape (D, 1)
bits = np.unpackbits(indices, axis=1, count=num_bits, bitorder='little').astype(float) # shape (D, num_bits)
# For each bit, expectation = sum_j [prob_j * bit_j]
y = probs.T @ bits # shape (T, num_bits)
# Tile the time points to have one x array per spin
x = np.tile(np.expand_dims(time_points, 1), (1, num_bits)) # shape (T, num_bits)
return x, y
def diagonalized_evolution(hamiltonian, initial_state, time, num_steps=100):
"""Diagonalize the given reduced Hamiltonian and evolve the initial state by exp(-i time*hamiltonian).
Args:
hamiltonian (np.ndarray(shape=(D, D), dtype=np.complex128)): Hamiltonian matrix divided by hbar.
initial_state (np.ndarray(shape=(D,), dtype=np.complex128)): Initial state vector.
time (float): Evolution time.
num_steps (int): Number of steps (T) to divide time into.
Returns:
np.ndarray(shape=(T,), dtype=float): Time points.
np.ndarray(shape=(D, T), dtype=np.complex128): State vector as a function of time.
"""
num_dim = hamiltonian.shape[0]
num_qubits = np.round(np.log2(num_dim)).astype(int)
# Create the array of time points
time_points = np.linspace(0., time, num_steps, endpoint=True)
## Diagonalize the Hamiltonian
eigvals, eigvectors = np.linalg.eigh(hamiltonian)
## Decompose the initial state vector into a linear combination of eigenvectors
# Matrix eigvectors has the form [v_0 v_1 v_2 ..], where v_i . v_j = delta_ij
# -> eigvectors^dagger @ initial_state = coefficients for the eigenvector decomposition of the initial state vector
initial_coeff = eigvectors.T.conjugate() @ initial_state
# Initial state as a matrix [c_0 v_0, c_1 v_1, ...] (shape (D, D))
initial_state_matrix = eigvectors * initial_coeff
## Time-evolve the initial state to each time point
# Phase at each time point (shape (D, T))
phase = np.outer(-1.j * eigvals, time_points)
phase_factor = np.exp(phase)
statevectors = initial_state_matrix @ phase_factor # shape (D, T)
return time_points, statevectors
def bit_expectations_counts(time_points, counts_list, num_bits):
"""Compute the bit expectation values from experiment results.
Args:
time_points (np.ndarray(shape=(T,), dtype=float)): Time points.
counts_list (List(Dict)): List (length T) of quantum experiment results, as given by Qiskit job.result().get_counts()
num_bits (int): Number of qubits
Returns:
np.ndarray(shape=(nstep, num_bits), dtype=float): Time points tiled for each bit.
np.ndarray(shape=(nstep, num_bits), dtype=float): Bit expectation values.
"""
if num_bits > 8:
raise NotImplementedError('Function not compatible with number of qubits > 8')
num_steps = len(counts_list)
x = np.tile(np.expand_dims(time_points, axis=1), (1, num_bits)) # shape (T, num_bits)
y = np.zeros_like(x)
for istep, counts in enumerate(counts_list):
counts = counts_list[istep]
total = 0
for bitstring, count in counts.items():
# 1. reverse the bitstring (last bit is the least significant)
# 2. map all bits to integers
# 3. convert to array
bits = np.array(list(map(int, reversed(bitstring))), dtype=float)
y[istep] += count * bits
total += count
y[istep] /= total
return x, y
def plot_heisenberg_spins(counts_list, num_spins, initial_state, omegadt, add_theory_curve=False, spin_component='z'):
"""Compute the expectation value of the Z(/X/Y) component of each spin in the Heisenberg model from the quantum
measurement results.
Args:
counts_list (List(Dict)): List of quantum experiment results, as given by Qiskit job.result().get_counts()
num_spins (int): Number of spins in the system.
initial_state (np.ndarray(shape=(2 ** num_spins), dtype=np.complex128)): Initial state vector.
omegadt (float): Hamiltonian parameter (H = -0.5 hbar omega sum_j [xx + yy + zz]) times time step.
add_theory_curve (bool): If True, compute the exact (non-Trotter) solution.
spin_component (str): Spin component to plot. Values 'x', 'y', or 'z'. Only affects the theory curve.
"""
# Number of steps
num_steps = len(counts_list)
# Figure and axes for the plot
fig, ax = plt.subplots(1, 1)
legend_items = []
legend_labels = []
if add_theory_curve:
# Construct the numerical Hamiltonian matrix from a list of Pauli operators
paulis = list()
for j in range(num_spins - 1):
paulis.append('I' * (num_spins - j - 2) + 'XX' + 'I' * j)
paulis.append('I' * (num_spins - j - 2) + 'YY' + 'I' * j)
paulis.append('I' * (num_spins - j - 2) + 'ZZ' + 'I' * j)
hamiltonian = SparsePauliOp(paulis).to_matrix()
# Compute the statevector as a function of time from Hamiltonian diagonalization
time_points, statevectors = diagonalized_evolution(-0.5 * hamiltonian, initial_state, omegadt * num_steps)
spin_basis_change = None
if spin_component == 'x':
spin_basis_change = np.array([[1., 1.], [1., -1.]], dtype=np.complex128) * np.sqrt(0.5)
elif spin_component == 'y':
spin_basis_change = np.array([[1., -1.j], [-1.j, 1.]], dtype=np.complex128) * np.sqrt(0.5)
if spin_basis_change is not None:
basis_change = tensor_product([spin_basis_change] * num_spins)
statevectors = basis_change @ statevectors
initial_state = basis_change @ initial_state
x, y = bit_expectations_sv(time_points, statevectors)
# Convert the bit expectations ([0, 1]) to spin expectations ([1, -1])
y = 1. - 2. * y
# Plot
lines = ax.plot(x, y)
colors = list(line.get_color() for line in lines)
dummy_line = mpl.lines.Line2D([0], [0])
dummy_line.update_from(lines[0])
dummy_line.set_color('black')
legend_items.append(dummy_line)
legend_labels.append('exact')
else:
colors = None
# Time points
time_points = np.linspace(0., num_steps * omegadt, num_steps + 1, endpoint=True)
# Prepend the initial "counts" to the experiment results
initial_probs = np.square(np.abs(initial_state))
fmt = f'{{:0{num_spins}b}}'
initial_counts = dict((fmt.format(idx), prob) for idx, prob in enumerate(initial_probs) if prob != 0.)
counts_list = [initial_counts] + counts_list
# Compute the bit expectation values from the counts
x, y = bit_expectations_counts(time_points, counts_list, num_spins)
# Convert the bit expectations ([0, 1]) to spin expectations ([1, -1])
y = 1. - 2. * y
# Plot
markers = ax.plot(x, y, 'o')
if colors is not None:
for marker, color in zip(markers, colors):
marker.set_color(color)
legend_items += markers
legend_labels += ['bit%d' % i for i in range(num_spins)]
ax.legend(legend_items, legend_labels)
ax.set_xlabel(r'$\omega t$')
ax.set_ylabel(r'$\langle S_z \rangle$')
def tensor_product(ops):
"""Recursively apply np.kron to construct a tensor product of the operators.
Args:
ops (List): List of (2, 2) arrays.
Returns:
np.ndarray(shape=(2 ** nops, 2 ** nops), dtype=np.complex128): Tensor product of ops.
"""
prod = 1.
for op in ops:
prod = np.kron(op, prod)
return prod