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

tobi update

parent 72b26dbb
import numpy as np
import math
import time
import matplotlib.pyplot as plt
class ObservationStation:
......@@ -11,20 +12,25 @@ class ObservationStation:
self.a = a
def execute(self, phi, r):
c = math.cos(2 * (self.a - phi))
s = math.sin(2 * (self.a - phi))
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():
phi1 = np.random.rand() * 2 * np.pi
phi2 = phi1 + np.pi / 2
return phi1, phi2
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 simulate():
def experiment():
pass
M = 100
def simulate(apply_time_window=False):
M = 1 # 100
N = int(1e6)
a = np.array([0, 1/4, 1/8, 3/8]) * np.pi
......@@ -36,33 +42,28 @@ def simulate():
station = (ObservationStation(T), ObservationStation(T))
# conducts EPRB experiment M times
for _ in range(M):
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])
# simulates 10^6 events
for n in range(N):
#if n%1000 == 0:
# print(p, n)
phi = generate_phis()
for s in range(2):
r = np.random.rand()
x, t_star = station[s].execute(phi[s], r)
results[n, p, s] = x
# 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
......@@ -70,13 +71,39 @@ def simulate():
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'|X| <= 2 is violated with |x| = {X_abs.max()}')
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'|X| <= 2 is satisfied with |x| = {X_abs.max()}')
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')
......
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