Skip to content

Commit 9828eeb

Browse files
mhjensenmhjensen
committed
adding codes
1 parent fe83636 commit 9828eeb

File tree

2 files changed

+337
-0
lines changed

2 files changed

+337
-0
lines changed

‎doc/src/week6/programs/onequbit.py‎

Lines changed: 143 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,143 @@
1+
import numpy as np
2+
import matplotlib.pyplot as plt
3+
import seaborn as sns; sns.set_theme(font_scale=1.5)
4+
from tqdm import tqdm
5+
6+
sigma_x = np.array([[0, 1], [1, 0]])
7+
sigma_y = np.array([[0, -1j], [1j, 0]])
8+
sigma_z = np.array([[1, 0], [0, -1]])
9+
I = np.eye(2)
10+
11+
def Hamiltonian(lmb):
12+
E1 = 0
13+
E2 = 4
14+
V11 = 3
15+
V22 = -3
16+
V12 = 0.2
17+
V21 = 0.2
18+
19+
eps = (E1 + E2) / 2
20+
omega = (E1 - E2) / 2
21+
c = (V11 + V22) / 2
22+
omega_z = (V11 - V22) / 2
23+
omega_x = V12
24+
25+
H0 = eps * I + omega * sigma_z
26+
H1 = c * I + omega_z * sigma_z + omega_x * sigma_x
27+
return H0 + lmb * H1
28+
29+
lmbvalues_ana = np.arange(0, 1, 0.01)
30+
eigvals_ana = np.zeros((len(lmbvalues_ana), 2))
31+
for index, lmb in enumerate(lmbvalues_ana):
32+
H = Hamiltonian(lmb)
33+
eigen, eigvecs = np.linalg.eig(H)
34+
permute = eigen.argsort()
35+
eigvals_ana[index] = eigen[permute]
36+
eigvecs = eigvecs[:,permute]
37+
38+
39+
fig, axs = plt.subplots(1, 1, figsize=(10, 10))
40+
for i in range(2):
41+
axs.plot(lmbvalues_ana, eigvals_ana[:,i], label=f'$E_{i+1}$')
42+
axs.set_xlabel(r'$\lambda$')
43+
axs.set_ylabel('Energy')
44+
axs.legend()
45+
plt.show()
46+
47+
#This was the standard eigenvalue problem. Let us now switch to our own implementation of the VQE.
48+
49+
from qc import *
50+
51+
def prepare_state(theta, phi, target = None):
52+
I = np.eye(2)
53+
sigma_x = np.array([[0, 1], [1, 0]])
54+
sigma_y = np.array([[0, -1j], [1j, 0]])
55+
state = np.array([1, 0])
56+
Rx = np.cos(theta/2) * I - 1j * np.sin(theta/2) * sigma_x
57+
Ry = np.cos(phi/2) * I - 1j * np.sin(phi/2) * sigma_y
58+
state = Ry @ Rx @ state
59+
if target is not None:
60+
state = target
61+
return state
62+
63+
def get_energy(angles, lmb, number_shots, target = None):
64+
theta, phi = angles[0], angles[1]
65+
# print(f'Theta = {theta}, Phi = {phi}')
66+
E1 = 0; E2 = 4; V11 = 3; V22 = -3; V12 = 0.2; V21 = 0.2
67+
68+
eps = (E1 + E2) / 2; omega = (E1 - E2) / 2; c = (V11 + V22) / 2; omega_z = (V11 - V22) / 2; omega_x = V12
69+
70+
init_state = prepare_state(theta, phi, target)
71+
qubit = One_qubit()
72+
qubit.set_state(init_state)
73+
measure_z = qubit.measure(number_shots)
74+
75+
qubit.set_state(init_state)
76+
qubit.apply_hadamard()
77+
measure_x = qubit.measure(number_shots)
78+
79+
# expected value of Z = (number of 0 measurements - number of 1 measurements)/ number of shots
80+
# number of 1 measurements = sum(measure_z)
81+
exp_val_z = (omega + lmb*omega_z)*(number_shots - 2*np.sum(measure_z)) / number_shots
82+
exp_val_x = lmb*omega_x*(number_shots - 2*np.sum(measure_x)) / number_shots
83+
exp_val_i = (eps + c*lmb)
84+
exp_val = (exp_val_z + exp_val_x + exp_val_i)
85+
return exp_val
86+
87+
def minimize_energy(lmb, number_shots, angles_0, learning_rate, max_epochs):
88+
# angles = np.random.uniform(low = 0, high = np.pi, size = 2)
89+
angles = angles_0 #lmb*np.array([np.pi, np.pi])
90+
epoch = 0
91+
delta_energy = 1
92+
energy = get_energy(angles, lmb, number_shots)
93+
while (epoch < max_epochs) and (delta_energy > 1e-4):
94+
grad = np.zeros_like(angles)
95+
for idx in range(angles.shape[0]):
96+
angles_temp = angles.copy()
97+
angles_temp[idx] += np.pi/2
98+
E_plus = get_energy(angles_temp, lmb, number_shots)
99+
angles_temp[idx] -= np.pi
100+
E_minus = get_energy(angles_temp, lmb, number_shots)
101+
grad[idx] = (E_plus - E_minus)/2
102+
angles -= learning_rate*grad
103+
new_energy = get_energy(angles, lmb, number_shots)
104+
delta_energy = np.abs(new_energy - energy)
105+
energy = new_energy
106+
epoch += 1
107+
return angles, epoch, (epoch < max_epochs), energy, delta_energy
108+
109+
110+
number_shots_search = 10_000
111+
number_shots = 10_000
112+
learning_rate = 0.3
113+
max_epochs = 400
114+
lmbvalues = np.linspace(0.0, 1.0, 30)
115+
min_energy = np.zeros(len(lmbvalues))
116+
epochs = np.zeros(len(lmbvalues))
117+
for index, lmb in enumerate(tqdm(lmbvalues)):
118+
memory = 0
119+
angles_0 = np.random.uniform(low = 0, high = np.pi, size = 2)
120+
angles, epochs[index], converged, energy, delta_energy = minimize_energy(lmb, number_shots_search, angles_0, learning_rate, max_epochs)
121+
if epochs[index] < (epochs[index-1] - 5):
122+
angles_0 = np.random.uniform(low = 0, high = np.pi, size = 2)
123+
angles, epochs[index], converged, energy, delta_energy = minimize_energy(lmb, number_shots_search, angles_0, learning_rate, max_epochs)
124+
min_energy[index] = get_energy(angles, lmb, number_shots)
125+
126+
from scipy.optimize import minimize
127+
number_shots = 10_000
128+
lmbvalues_scipy = np.linspace(0.0, 1.0, 50)
129+
min_energy_scipy = np.zeros(len(lmbvalues_scipy))
130+
for index, lmb in enumerate(tqdm(lmbvalues_scipy)):
131+
angles_start = np.random.uniform(low = 0, high = np.pi, size = 4)
132+
res = minimize(get_energy, angles_start, args = (lmb, number_shots), method = 'Powell', options = {'maxiter': 1000}, tol = 1e-5)
133+
min_energy_scipy[index] = res.fun
134+
135+
fig, axs = plt.subplots(1, 1, figsize=(10, 10))
136+
for i in range(2):
137+
axs.plot(lmbvalues_ana, eigvals_ana[:,i], label=f'$E_{i+1}$', color = '#4c72b0')
138+
axs.scatter(lmbvalues, min_energy, label = 'VQE eigenvalues', color = '#dd8452')
139+
axs.scatter(lmbvalues_scipy, min_energy_scipy, label = 'VQE Scipy', color = '#55a868')
140+
axs.set_xlabel(r'$\lambda$')
141+
axs.set_ylabel('Energy')
142+
plt.legend()
143+
plt.show()

‎doc/src/week6/programs/qc.py‎

Lines changed: 194 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,194 @@
1+
import numpy as np
2+
3+
class One_qubit:
4+
def __init__(self):
5+
self.state = np.zeros(2, dtype=np.complex_)
6+
self.I = np.eye(2)
7+
self.Z = np.array([[1, 0], [0, -1]])
8+
self.X = np.array([[0, 1], [1, 0]])
9+
self.Y = np.array([[0, -1j], [1j, 0]])
10+
self.H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
11+
self.S = np.array([[1, 0], [0, 1j]])
12+
13+
def set_state(self, state):
14+
if abs(np.linalg.norm(state) - 1) > 1e-10:
15+
raise ValueError("The state vector must be normalized.")
16+
self.state = state
17+
18+
def apply_hadamard(self):
19+
self.state = np.dot(self.H, self.state)
20+
return self.state
21+
22+
def apply_x(self):
23+
self.state = np.dot(self.X, self.state)
24+
return self.state
25+
26+
def apply_y(self):
27+
self.state = np.dot(self.Y, self.state)
28+
return self.state
29+
30+
def apply_z(self):
31+
self.state = np.dot(self.Z, self.state)
32+
return self.state
33+
34+
def measure(self, num_shots=1):
35+
prob = np.abs(self.state)**2
36+
possible = np.arange(len(self.state)) #possible measurement outcomes
37+
outcome = np.random.choice(possible, p=prob, size = num_shots) #measurement outcome
38+
self.state = np.zeros_like(self.state) #set state to the measurement outcome
39+
self.state[outcome[-1]] = 1
40+
return outcome
41+
42+
def rotate_x(self, theta):
43+
# implement rotation around x axis
44+
Rx = np.cos(theta/2) * self.I - 1j * np.sin(theta/2) * self.X
45+
self.state = np.dot(Rx, self.state)
46+
47+
def rotate_y(self, phi):
48+
# implement rotation around y axis
49+
Ry = np.cos(phi/2) * self.I - 1j * np.sin(phi/2) * self.Y
50+
self.state = np.dot(Ry, self.state)
51+
52+
class Two_qubit(One_qubit):
53+
def __init__(self):
54+
super().__init__()
55+
self.state = np.zeros(4, dtype=np.complex_)
56+
self.CNOT01 = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])
57+
self.CNOT10 = np.array([[1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0]])
58+
self.SWAP = np.array([[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]])
59+
60+
# np.random.seed(0)
61+
62+
def apply_cnot01(self):
63+
self.state = np.dot(self.CNOT01, self.state)
64+
return self.state
65+
66+
def apply_cnot10(self):
67+
self.state = np.dot(self.CNOT10, self.state)
68+
return self.state
69+
70+
def apply_swap(self):
71+
self.state = np.dot(self.SWAP, self.state)
72+
return self.state
73+
74+
def apply_hadamard(self, qubit):
75+
if qubit == 0:
76+
self.state = np.kron(self.H, self.I).dot(self.state)
77+
elif qubit == 1:
78+
self.state = np.kron(self.I, self.H).dot(self.state)
79+
return self.state
80+
81+
def apply_sdag(self, qubit):
82+
if qubit == 0:
83+
self.state = np.kron(self.S.conj().T, self.I).dot(self.state)
84+
elif qubit == 1:
85+
self.state = np.kron(self.I, self.S.conj().T).dot(self.state)
86+
return self.state
87+
88+
def apply_x(self, qubit):
89+
if qubit == 0:
90+
self.state = np.kron(self.X, self.I).dot(self.state)
91+
elif qubit == 1:
92+
self.state = np.kron(self.I, self.X).dot(self.state)
93+
return self.state
94+
95+
def apply_y(self, qubit):
96+
if qubit == 0:
97+
self.state = np.kron(self.Y, self.I).dot(self.state)
98+
elif qubit == 1:
99+
self.state = np.kron(self.I, self.Y).dot(self.state)
100+
return self.state
101+
102+
def apply_z(self, qubit):
103+
if qubit == 0:
104+
self.state = np.kron(self.Z, self.I).dot(self.state)
105+
elif qubit == 1:
106+
self.state = np.kron(self.I, self.Z).dot(self.state)
107+
return self.state
108+
109+
def rotate_x(self, theta, qubit):
110+
Rx = np.cos(theta/2) * self.I - 1j * np.sin(theta/2) * self.X
111+
if qubit == 0:
112+
self.state = np.kron(Rx, self.I).dot(self.state)
113+
elif qubit == 1:
114+
self.state = np.kron(self.I, Rx).dot(self.state)
115+
return self.state
116+
117+
def rotate_y(self, phi, qubit):
118+
# implement rotation around y axis
119+
Ry = np.cos(phi/2) * self.I - 1j * np.sin(phi/2) * self.Y
120+
if qubit == 0:
121+
self.state = np.kron(Ry, self.I).dot(self.state)
122+
elif qubit == 1:
123+
self.state = np.kron(self.I, Ry).dot(self.state)
124+
return self.state
125+
126+
class Four_qubit(Two_qubit):
127+
#This is based on the two qubit class since the Lipkin model at most acts on two qubits at the time
128+
def __init__(self):
129+
super().__init__()
130+
self.state = np.zeros(16, dtype=np.complex_)
131+
132+
def apply_cnot10(self, qubit1):
133+
# can only be applied for adjecent qubits
134+
if qubit1 == 0:
135+
op = np.kron(self.CNOT10, np.kron(self.I, self.I))
136+
return np.dot(op, self.state)
137+
elif qubit1 == 1:
138+
op = np.kron(self.I, np.kron(self.CNOT10, self.I))
139+
return np.dot(op, self.state)
140+
elif qubit1 == 2:
141+
op = np.kron(self.I, np.kron(self.I, self.CNOT10))
142+
return np.dot(op, self.state)
143+
else:
144+
print('qubit1 must be 0, 1, or 2')
145+
146+
def apply_cnot01(self, qubit1):
147+
# can only be applied for adjecent qubits
148+
if qubit1 == 0:
149+
op = np.kron(self.CNOT01, np.kron(self.I, self.I))
150+
return np.dot(op, self.state)
151+
elif qubit1 == 1:
152+
op = np.kron(self.I, np.kron(self.CNOT01, self.I))
153+
return np.dot(op, self.state)
154+
elif qubit1 == 2:
155+
op = np.kron(self.I, np.kron(self.I, self.CNOT01))
156+
return np.dot(op, self.state)
157+
else:
158+
print('qubit1 must be 0, 1, or 2')
159+
160+
def apply_swap(self, qubit1):
161+
# can only be applied for adjecent qubits
162+
if qubit1 == 0:
163+
op = np.kron(self.SWAP, np.kron(self.I, self.I))
164+
return np.dot(op, self.state)
165+
elif qubit1 == 1:
166+
op = np.kron(self.I, np.kron(self.SWAP, self.I))
167+
return np.dot(op, self.state)
168+
elif qubit1 == 2:
169+
op = np.kron(self.I, np.kron(self.I, self.SWAP))
170+
return np.dot(op, self.state)
171+
else:
172+
print('qubit1 must be 0, 1, or 2')
173+
174+
def apply_hadamard(self, qubit):
175+
if qubit == 0:
176+
self.state = np.kron(self.H, np.kron(self.I, np.kron(self.I, self.I))).dot(self.state)
177+
elif qubit == 1:
178+
self.state = np.kron(self.I, np.kron(self.H, np.kron(self.I, self.I))).dot(self.state)
179+
elif qubit == 2:
180+
self.state = np.kron(self.I, np.kron(self.I, np.kron(self.H, self.I))).dot(self.state)
181+
elif qubit == 3:
182+
self.state = np.kron(self.I, np.kron(self.I, np.kron(self.I, self.H))).dot(self.state)
183+
return self.state
184+
185+
def apply_sdag(self, qubit):
186+
if qubit == 0:
187+
self.state = np.kron(self.S.conj().T, np.kron(self.I, np.kron(self.I, self.I))).dot(self.state)
188+
elif qubit == 1:
189+
self.state = np.kron(self.I, np.kron(self.S.conj().T, np.kron(self.I, self.I))).dot(self.state)
190+
elif qubit == 2:
191+
self.state = np.kron(self.I, np.kron(self.I, np.kron(self.S.conj().T, self.I))).dot(self.state)
192+
elif qubit == 3:
193+
self.state = np.kron(self.I, np.kron(self.I, np.kron(self.I, self.S.conj().T))).dot(self.state)
194+
return self.state

0 commit comments

Comments
 (0)