Commit bbcaf66f authored by Tobias Jörg Theodor Kempe's avatar Tobias Jörg Theodor Kempe
Browse files

resolves merge conflicts

parents 451e259b dfa6553c
import numpy as np
<<<<<<< HEAD
import time
import matplotlib.pyplot as plt
class ObservationStation:
def __init__(self, T):
self.T = T
def set_a(self, a):
self.a = a
def execute(self, phi, r):
c = np.cos(2 * (self.a - phi))
s = np.sin(2 * (self.a - phi))
x = np.sign(c)
t_star = r * self.T * s ** 2
return x, t_star
def generate_phis(count = 1):
phis = np.empty((count, 4, 2), dtype=np.float64)
phis[:, :, 0] = np.random.rand(count, 4) * 2 * np.pi
phis[:, :, 1] = phis[:, :, 0] + np.pi / 2
return phis
def experiment():
pass
def simulate(apply_time_window=False):
M = 1 # 100
N = int(1e6)
a = np.array([0, 1/4, 1/8, 3/8]) * np.pi
W = 1
T = 1000
station = (ObservationStation(T), ObservationStation(T))
# conducts EPRB experiment M times
for m in range(M):
results = np.zeros((N, 4, 2), dtype=np.int8)
phis = generate_phis(N)
rs = np.random.rand(N, 4, 2)
# iterates through all pairs of settings
for p in range(4):
# selects a unique pair of settings from all possible settings
station[0].set_a(a[int(p/2)])
station[1].set_a(a[2+p%2])
# alternative approach to simulating 10^6 events
x1, t_star1 = station[0].execute(phis[:, p, 0], rs[:, p, 0])
results[:, p, 0] = x1
x2, t_star2 = station[1].execute(phis[:, p, 1], rs[:, p, 1])
results[:, p, 1] = x2
# time window selection
# calculates absolute X
products = np.prod(results, axis=2)
products[:, 2] *= -1
sums = np.sum(products, axis=1)
X_abs = np.abs(sums)
# calculates expectancy values of x(a1) and x(a2)
temp = np.mean(results, axis=0)
E_1 = temp[:, 0]
E_2 = temp[:, 1]
E_1_times_E_2 = E_1 * E_2
# calculates expectancy value of x(a1)*x(a2)
temp = np.prod(results, axis=2)
E_12 = np.mean(temp, axis=0)
# plots the results
# do we plot the results for each configuration separately?
plt.figure()
plt.plot(, E_1_times_E_2p[0])
if np.any(X_abs >= 2):
print(f'{m} - |X| <= 2 is violated with |x| = {X_abs.max()}')
else:
print(f'{m} - |X| <= 2 is satisfied with |x| = {X_abs.max()}')
if np.mean(X_abs) >= 2:
print(f'{m} - |X| <= 2 is violated with |x| = {np.mean(X_abs)}')
else:
print(f'{m} - |X| <= 2 is satisfied with |x| = {np.mean(X_abs)}')
start = time.time()
simulate()
duration = time.time() - start
print(f'total simulation duration: {round(duration, 2)} s')
=======
import matplotlib.pyplot as plt
import time
# settings
N = int(100000)
T_0 = 1000
W = 1
showPlots = True
def simulate(phi1, phi2):
# generate random phis
thetas = np.empty((N, 2), dtype=np.float32)
thetas[:, 0] = np.random.rand(N) * 2 * np.pi
thetas[:, 1] = thetas[:, 0] + np.pi / 2
# calculate outcome of observation station
x1 = (np.cos(2 * (phi1 - thetas[:, 0])) >= (2*np.random.rand(N) -1)).reshape(N, 1)
t1 = T_0 * np.random.rand(N) * np.sin(2 * (phi1 - thetas[:, 0])) ** 4
x2 = (np.cos(2 * (phi2 - thetas[:, 1])) >= (2*np.random.rand(N) -1)).reshape(N, 1)
t2 = T_0 * np.random.rand(N) * np.sin(2 * (phi2 - thetas[:, 1])) ** 4
# determine which events are conincidences (1) and which not (0)
coincidences = np.abs(t1 - t2) < W
# zip measurement outcome together
x1x2 = np.concatenate((x1, x2), axis=1)
counts = np.empty((2, 4), dtype=np.int32)
counts[0, 0] = np.sum(np.all(x1x2 == np.array([0, 0]), axis=1), axis=0)
counts[0, 1] = np.sum(np.all(x1x2 == np.array([0, 1]), axis=1), axis=0)
counts[0, 2] = np.sum(np.all(x1x2 == np.array([1, 0]), axis=1), axis=0)
counts[0, 3] = np.sum(np.all(x1x2 == np.array([1, 1]), axis=1), axis=0)
counts[1, 0] = np.sum(np.logical_and(np.all(x1x2 == np.array([0, 0]), axis=1), coincidences), axis=0)
counts[1, 1] = np.sum(np.logical_and(np.all(x1x2 == np.array([0, 1]), axis=1), coincidences), axis=0)
counts[1, 2] = np.sum(np.logical_and(np.all(x1x2 == np.array([1, 0]), axis=1), coincidences), axis=0)
counts[1, 3] = np.sum(np.logical_and(np.all(x1x2 == np.array([1, 1]), axis=1), coincidences), axis=0)
E_all = np.sum(counts, axis=1).astype(float)
E_12 = (counts[:, 0] - counts[:, 1] - counts[:, 2] + counts[:, 3]).astype(float)
E_1 = (counts[:, 0] + counts[:, 1] - counts[:, 2] - counts[:, 3]).astype(float)
E_2 = (counts[:, 0] - counts[:, 1] + counts[:, 2] - counts[:, 3]).astype(float)
if (np.all(E_all > 0)):
E_12 /= E_all
E_1 /= E_all
E_2 /= E_all
return E_12, E_1, E_2
def try_all_configs(n_configs):
E_12 = np.empty((n_configs, 2), dtype=np.float32)
E_1 = np.empty((n_configs, 2), dtype=np.float32)
E_2 = np.empty((n_configs, 2), dtype=np.float32)
for i in range(n_configs):
phi1 = 0
phi2 = i * 2 * np.pi / n_configs
E_12_, E_1_, E_2_ = simulate(phi1, phi2)
E_12[i, :] = E_12_
E_1[i, :] = E_1_
E_2[i, :] = E_2_
# correlation
phi = np.linspace(0, 360/n_configs*(n_configs-1), n_configs)
phi_theory = np.linspace(0, 360/n_configs*(n_configs-1), n_configs*100)
theory = -np.cos(2 * phi_theory * (2 * np.pi / 360))
plt.figure()
plt.xlabel(r'$\Delta \varphi$')
plt.ylabel(r'$<x_1 x_2>$')
plt.plot(phi_theory, theory, '.', markersize=1, color='orange')
plt.plot(phi, E_12[:, 0], 'o', color='blue')
#plt.savefig('correlation.pdf')
if showPlots:
plt.show()
plt.close()
# product of expectation values
line = np.zeros(phi_theory.shape)
plt.figure()
plt.xlabel(r'$\Delta \varphi$')
plt.ylabel(r'$<x_1><x_2>$')
plt.ylim([-1, 1])
plt.xticks([0, 45, 90, 135, 180, 225, 270, 305, 360])
plt.plot(phi_theory, line, '.', markersize=1, color='orange')
plt.plot(phi, E_1[:, 0]*E_2[:, 0], 'o', color='blue')
#plt.savefig('expectation_value.pdf')
if showPlots:
plt.show()
plt.close()
# correlation
phi = np.linspace(0, 360/n_configs*(n_configs-1), n_configs)
phi_theory = np.linspace(0, 360/n_configs*(n_configs-1), n_configs*100)
theory = -np.cos(2 * phi_theory * (2 * np.pi / 360))
plt.figure()
plt.xlabel(r'$\Delta \varphi$')
plt.ylabel(r'$<x_1 x_2>$')
plt.plot(phi_theory, theory, '.', markersize=1, color='orange')
plt.plot(phi, E_12[:, 1], 'o', color='blue')
#plt.savefig('correlation.pdf')
if showPlots:
plt.show()
plt.close()
# product of expectation values
line = np.zeros(phi_theory.shape)
plt.figure()
plt.xlabel(r'$\Delta \varphi$')
plt.ylabel(r'$<x_1><x_2>$')
plt.ylim([-1, 1])
plt.xticks([0, 45, 90, 135, 180, 225, 270, 305, 360])
plt.plot(phi_theory, line, '.', markersize=1, color='orange')
plt.plot(phi, E_1[:, 1]*E_2[:, 1], 'o', color='blue')
#plt.savefig('expectation_value.pdf')
if showPlots:
plt.show()
plt.close()
start = time.time()
try_all_configs(32)
duration = time.time() - start
print(f'simulation took {round(duration, 2)} seconds')
>>>>>>> 306f34e3cfd8a671917195d947d9ef22046bb21a
import numpy as np
import matplotlib.pyplot as plt
# some initialization
wxyz = 1234
np.random.seed(wxyz)
# input data
nsteps = 32
nsamples = 100000
T0 = 1000 # (ns), maximum time delay
W = 1 # (ns), time coincidence window
HWP2 = 0 # degrees, orientation of wave plate 2 (EOM2)
cHWP2 = np.cos(HWP2*np.pi/180)
sHWP2 = np.sin(HWP2*np.pi/180)
def analyzer(c, s, cHWP, sHWP, T0):
# EOM : plane rotation
c2 = cHWP*c + sHWP*s
s2 = -sHWP*c + cHWP*s
x = c2*c2 - s2*s2 # cos(2(x-a))
y = 2*c2*s2 # sin(2(x-a))
# Malus law
r0 = np.random.rand()
if(x>2*r0 - 1):
j = 0 # +1 event
else:
j = 1 #-1 event
# time delay
r1 = np.random.rand()
l = y**4*T0*r1 # delay time: T_0 sin(2*(theta1 - x))**4
return j, l
count = np.zeros(shape=(2, 2, 2, nsteps), dtype=np.int32) # set all counts to zero
for ipsi0 in range(nsteps): # loop over different settings of EOM1
cHWP1 = np.cos(ipsi0*2*np.pi/nsteps)
sHWP1 = np.sin(ipsi0*2*np.pi/nsteps)
for i in range(nsamples):
# source
r0 = np.random.rand()
c1 = np.cos(r0*2*np.pi) # polarization angle x of particle going to station 1
s1 = np.sin(r0*2*np.pi) # polarization angle x + pi/2 of particle going to station 2
c2 = -s1
s2 = c1
# station 1
j1, l1 = analyzer(c1, s1, cHWP1, sHWP1, T0)
# station 2
j2, l2 = analyzer(c2, s2, cHWP2, sHWP2, T0)
# count
count[j1, j2, 0, ipsi0] = count[j1, j2, 0, ipsi0] + 1 # Malus law model
if(abs(l1 - l2)<W):
count[j1, j2, 1, ipsi0] = count[j1, j2, 1, ipsi0] + 1 # Malus law model + time window
# data analysis
tot = np.zeros(shape=(2, nsteps), dtype=np.int32)
E12 = np.zeros(shape=(2, nsteps), dtype=np.float64)
E1 = np.zeros(shape=(2, nsteps), dtype=np.float64)
E2 = np.zeros(shape=(2, nsteps), dtype=np.float64)
for j in range(nsteps):
for i in [0, 1]:
# i = 0 <==> no time window, i = 1 <=> use time coincidences
#tot[i, j] = np.sum(count[:, :, i, j])
tot[i, j] = count[0, 0, i, j] + count[1, 1, i, j] + count[1, 0, i, j] + count[0, 1, i, j]
E12[i, j] = count[0, 0, i, j] + count[1, 1, i, j] - count[1, 0, i, j] - count[0, 1, i, j]
E1[i, j] = count[0, 0, i, j] + count[0, 1, i, j] - count[1, 1, i, j] - count[1, 0, i, j]
E2[i, j] = count[0, 0, i, j] + count[1, 0, i, j] - count[1, 1, i, j] - count[0, 1, i, j]
if(tot[i, j]>0):
E12[i, j] = E12[i, j]/tot[i, j]
E1[i, j] = E1[i, j]/tot[i, j]
E2[i, j] = E2[i, j]/tot[i, j]
theta = np.linspace(0, 360, nsteps)
theta_theory = np.linspace(0, 360, nsteps*100)
theory = []
for j in range(nsteps*100):
theory.append(-np.cos(2 * j * 2 * np.pi / nsteps / 100))
plt.plot(theta, E12[0, :], 'o')
plt.plot(theta_theory, theory, '.', markersize=1)
plt.savefig('e12.pdf')
plt.close()
\ No newline at end of file
import numpy as np
<<<<<<< HEAD
import time
import matplotlib.pyplot as plt
class ObservationStation:
def __init__(self, T):
self.T = T
def set_a(self, a):
self.a = a
def execute(self, phi, r):
c = np.cos(2 * (self.a - phi))
s = np.sin(2 * (self.a - phi))
x = np.sign(c)
t_star = r * self.T * s ** 2
return x, t_star
def generate_phis(count = 1):
phis = np.empty((count, 4, 2), dtype=np.float64)
phis[:, :, 0] = np.random.rand(count, 4) * 2 * np.pi
phis[:, :, 1] = phis[:, :, 0] + np.pi / 2
return phis
def experiment():
pass
def simulate(apply_time_window=False):
M = 1 # 100
N = int(1e6)
a = np.array([0, 1/4, 1/8, 3/8]) * np.pi
W = 1
T = 1000
station = (ObservationStation(T), ObservationStation(T))
# conducts EPRB experiment M times
for m in range(M):
results = np.zeros((N, 4, 2), dtype=np.int8)
phis = generate_phis(N)
rs = np.random.rand(N, 4, 2)
# iterates through all pairs of settings
for p in range(4):
# selects a unique pair of settings from all possible settings
station[0].set_a(a[int(p/2)])
station[1].set_a(a[2+p%2])
# alternative approach to simulating 10^6 events
x1, t_star1 = station[0].execute(phis[:, p, 0], rs[:, p, 0])
results[:, p, 0] = x1
x2, t_star2 = station[1].execute(phis[:, p, 1], rs[:, p, 1])
results[:, p, 1] = x2
# time window selection
# calculates absolute X
products = np.prod(results, axis=2)
products[:, 2] *= -1
sums = np.sum(products, axis=1)
X_abs = np.abs(sums)
# calculates expectancy values of x(a1) and x(a2)
temp = np.mean(results, axis=0)
E_1 = temp[:, 0]
E_2 = temp[:, 1]
E_1_times_E_2 = E_1 * E_2
# calculates expectancy value of x(a1)*x(a2)
temp = np.prod(results, axis=2)
E_12 = np.mean(temp, axis=0)
# plots the results
# do we plot the results for each configuration separately?
plt.figure()
plt.plot(, E_1_times_E_2p[0])
if np.any(X_abs >= 2):
print(f'{m} - |X| <= 2 is violated with |x| = {X_abs.max()}')
else:
print(f'{m} - |X| <= 2 is satisfied with |x| = {X_abs.max()}')
if np.mean(X_abs) >= 2:
print(f'{m} - |X| <= 2 is violated with |x| = {np.mean(X_abs)}')
else:
print(f'{m} - |X| <= 2 is satisfied with |x| = {np.mean(X_abs)}')
start = time.time()
simulate()
duration = time.time() - start
print(f'total simulation duration: {round(duration, 2)} s')
=======
import matplotlib.pyplot as plt
import time
......@@ -127,4 +237,5 @@ def try_all_configs(n_configs):
start = time.time()
try_all_configs(32)
duration = time.time() - start
print(f'simulation took {round(duration, 2)} seconds')
\ No newline at end of file
print(f'simulation took {round(duration, 2)} seconds')
>>>>>>> 306f34e3cfd8a671917195d947d9ef22046bb21a
<<<<<<< HEAD:python/original.py
import numpy as np
import matplotlib.pyplot as plt
import time
......@@ -164,4 +165,8 @@ plt.close()
# displays stopwatch
duration = time.time() - start
print(f'simulation took {round(duration, 2)} seconds')
\ No newline at end of file
print(f'simulation took {round(duration, 2)} seconds')
=======
import example
import eprb
>>>>>>> master:python/main.py
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment