optimized.py 4.55 KB
Newer Older
Tobias Jörg Theodor Kempe's avatar
2nd try    
Tobias Jörg Theodor Kempe committed
1
2
3
4
5
6

import numpy as np
import matplotlib.pyplot as plt
import time

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

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
20
    x1 = (np.cos(2 * (phi1 - thetas[:, 0])) >= (2*np.random.rand(N) -1)).reshape(N, 1)
21
    t1 = T_0 * np.random.rand(N) * np.sin(2 * (phi1 - thetas[:, 0])) ** 4
22
    x2 = (np.cos(2 * (phi2 - thetas[:, 1])) >= (2*np.random.rand(N) -1)).reshape(N, 1)
23
    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
24
25
26
27
28
29
30
31
32
33
34
35
36
37
                             
    # 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)
    
38
39
40
41
    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
42
43
44
45
46
47
48
49
50
51
52
53
54
55
    
    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
        
        
56
def try_all_configs(n_configs):
Tobias Jörg Theodor Kempe's avatar
2nd try    
Tobias Jörg Theodor Kempe committed
57
    
58
59
60
    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
61
    
62
    for i in range(n_configs):
Tobias Jörg Theodor Kempe's avatar
2nd try    
Tobias Jörg Theodor Kempe committed
63
        phi1 = 0
64
        phi2 = i * 2 * np.pi / n_configs
Tobias Jörg Theodor Kempe's avatar
2nd try    
Tobias Jörg Theodor Kempe committed
65
66
67
68
69
70
        
        E_12_, E_1_, E_2_ = simulate(phi1, phi2)
        E_12[i, :] = E_12_
        E_1[i, :] = E_1_
        E_2[i, :] = E_2_

71
    # correlation
72
73
    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)
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
    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()

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
    # 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()

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