optimized.py 5.29 KB
Newer Older
Tobias Jörg Theodor Kempe's avatar
2nd try    
Tobias Jörg Theodor Kempe committed
1
2
3
4
import numpy as np
import matplotlib.pyplot as plt
import time

Tobias Jörg Theodor Kempe's avatar
Tobias Jörg Theodor Kempe committed
5
6
7
8
# initialization
wxyz = 4754
np.random.seed(wxyz)

Tobias Jörg Theodor Kempe's avatar
2nd try    
Tobias Jörg Theodor Kempe committed
9
# settings
10
N = int(100000)
Tobias Jörg Theodor Kempe's avatar
2nd try    
Tobias Jörg Theodor Kempe committed
11
12
T_0 = 1000
W = 1
13
showPlots = True
Tobias Jörg Theodor Kempe's avatar
2nd try    
Tobias Jörg Theodor Kempe committed
14
15
16
17
18
19
20
21
22

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
23
    x1 = (np.cos(2 * (phi1 - thetas[:, 0])) >= (2*np.random.rand(N) -1)).reshape(N, 1)
24
    t1 = T_0 * np.random.rand(N) * np.sin(2 * (phi1 - thetas[:, 0])) ** 4
25
    x2 = (np.cos(2 * (phi2 - thetas[:, 1])) >= (2*np.random.rand(N) -1)).reshape(N, 1)
26
    t2 = T_0 * np.random.rand(N) * np.sin(2 * (phi2 - thetas[:, 1])) ** 4
Tobias Jörg Theodor Kempe's avatar
2nd try    
Tobias Jörg Theodor Kempe committed
27
28
29
30
31
32
33
                             
    # 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)
    
34
35
36
    # creates array to store counts in
    # 1st index ~ without time constraints (0) or with coincidence constraint (1)
    # 2nd index ~ the respective outcome [0, 0] (0), [0, 1] (1), [1, 0] (2), or [1, 1] (3)
Tobias Jörg Theodor Kempe's avatar
2nd try    
Tobias Jörg Theodor Kempe committed
37
38
    counts = np.empty((2, 4), dtype=np.int32)
    
39
    # calculates counts without time constraints
Tobias Jörg Theodor Kempe's avatar
2nd try    
Tobias Jörg Theodor Kempe committed
40
41
42
43
44
    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)
    
45
    # calculates counts with coincidence constraint
46
47
48
49
    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)
Tobias Jörg Theodor Kempe's avatar
2nd try    
Tobias Jörg Theodor Kempe committed
50
    
51
    # computes expectancy values for x1*x2, x1 and x2
Tobias Jörg Theodor Kempe's avatar
2nd try    
Tobias Jörg Theodor Kempe committed
52
    E_all = np.sum(counts, axis=1).astype(float)
53
    print(E_all)
Tobias Jörg Theodor Kempe's avatar
2nd try    
Tobias Jörg Theodor Kempe committed
54
55
56
57
58
59
60
61
62
63
64
    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
        
        
65
def try_all_configs(n_configs):
Tobias Jörg Theodor Kempe's avatar
2nd try    
Tobias Jörg Theodor Kempe committed
66
    
67
    # allocates storage for each required expectancy value
68
69
70
    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)
Tobias Jörg Theodor Kempe's avatar
2nd try    
Tobias Jörg Theodor Kempe committed
71
    
72
    # simulates all configurations
73
    for i in range(n_configs):
Tobias Jörg Theodor Kempe's avatar
2nd try    
Tobias Jörg Theodor Kempe committed
74
        phi1 = 0
75
        phi2 = i * 2 * np.pi / n_configs
Tobias Jörg Theodor Kempe's avatar
2nd try    
Tobias Jörg Theodor Kempe committed
76
77
78
79
80
81
        
        E_12_, E_1_, E_2_ = simulate(phi1, phi2)
        E_12[i, :] = E_12_
        E_1[i, :] = E_1_
        E_2[i, :] = E_2_

82
    # correlation
83
84
    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)
85
86
    theory = -np.cos(2 * phi_theory * (2 * np.pi / 360))
    plt.figure()
Orkun Şensebat's avatar
Orkun Şensebat committed
87
    plt.xlabel(r'$\Delta \varphi$ in degrees')
88
    plt.ylabel(r'$\left<x_1 x_2\right>$')
89
    plt.xticks([0, 45, 90, 135, 180, 225, 270, 305, 360])
90
91
    plt.plot(phi_theory, theory, '.', markersize=1, color='orange')
    plt.plot(phi, E_12[:, 0], 'o', color='blue')
92
    plt.savefig('correlation.pdf')
93
94
95
96
97
98
99
    if showPlots:
        plt.show()
    plt.close()
    
    # product of expectation values
    line = np.zeros(phi_theory.shape)
    plt.figure()
Orkun Şensebat's avatar
Orkun Şensebat committed
100
    plt.xlabel(r'$\Delta \varphi$ in degrees')
101
    plt.ylabel(r'$\left<x_1\right>\left<x_2\right>$')
102
103
104
105
    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')
106
    plt.savefig('expectation_value.pdf')
107
108
109
110
    if showPlots:
        plt.show()
    plt.close()

111
112
113
114
115
    # 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()
Orkun Şensebat's avatar
Orkun Şensebat committed
116
    plt.xlabel(r'$\Delta \varphi$ in degrees')
117
    plt.ylabel(r'$\left<x_1 x_2\right>$')
118
    plt.xticks([0, 45, 90, 135, 180, 225, 270, 305, 360])
119
120
    plt.plot(phi_theory, theory, '.', markersize=1, color='orange')
    plt.plot(phi, E_12[:, 1], 'o', color='blue')
121
    plt.savefig('correlation_coinc.pdf')
122
123
124
125
126
127
128
    if showPlots:
        plt.show()
    plt.close()
    
    # product of expectation values
    line = np.zeros(phi_theory.shape)
    plt.figure()
Orkun Şensebat's avatar
Orkun Şensebat committed
129
    plt.xlabel(r'$\Delta \varphi$ in degrees')
130
    plt.ylabel(r'$\left<x_1\right>\left<x_2\right>$')
131
132
133
134
    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')
135
    plt.savefig('expectation_value_coinc.pdf')
136
137
138
139
    if showPlots:
        plt.show()
    plt.close()

Tobias Jörg Theodor Kempe's avatar
2nd try    
Tobias Jörg Theodor Kempe committed
140
141
142
start = time.time()
try_all_configs(32)
duration = time.time() - start
143
print(f'simulation took {round(duration, 2)} seconds')