Simulated optics experiment/Data analysis
Simulated optics experiment/Data analysis is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
Introduction
In this task, you will write a program to analyze the "raw data" outputted by one of the Simulated optics experiment/Simulator implementations. The analysis will estimate the following numbers:
- Four correlation coefficients: one for each of the four possible combinations of the two polarizing beam splitter angles. These measure the interrelatedness of light detection events on the two sides of the experiment.
- A CHSH contrast, which is derived from the correlation coefficients.
The estimated CHSH contrast is what we are really interested in, because what we intend to show is that it contradicts what many physicists believe.
Quantum mechanics predicts a CHSH contrast of approximately . The Simulated optics experiment/Simulator simulation will be shown (by our analysis) to have a CHSH contrast very close to that. Note that our simulation does not employ any quantum mechanics, but instead simulates light "classically". Also our simulation clearly involves only causation that in each instance takes a finite amount of forward-moving time. That is, there is no "instantaneous action at a distance".
The prevailing view in physics, however, is that a CHSH contrast greater than 2 (which we certainly achieve) is possible only with quantum mechanics, and also that taking measurements on one side of the apparatus "instantaneously" determines what the measurements will be on the other side. Although quantum mechanics actually has no concept of "causation", this "determination" is often described as being "instantaneous action at a distance".
Thus, we are left with a conundrum: we clearly get a CHSH contrast not only greater than 2, but in fact very close to the theoretical value that quantum mechanics predicts. This is, in fact, a validation of the quantum mechanical mathematics: the mathematics is producing the correct result for this kind of experiment. On the other hand, the strong prevailing view in physics is that this cannot be happening if the physics involved is "classical"--which our simulated physics very much is. The prevailing view is that experiments similar to ours prove that "instantaneous action at a distance" is a real thing that exists. One sometimes sees such "proof" touted in the popular media.
I leave it to each individual to think about that conundrum and come to their own opinions. I am not a physicist and probably neither are you, so we can think whatever we want, and no one cares. But we likely are both computer programmers. What we mainly want to do here is demonstrate the use of computers to simulate physical events and analyze the data generated.
Task description
You will write a program to perform the mathematical data analysis described in the open access paper
A. F. Kracklauer, ‘EPR-B correlations: non-locality or geometry?’, J. Nonlinear Math. Phys. 11 (Supp.) 104–109 (2004). https://doi.org/10.2991/jnmp.2004.11.s1.13
The program will read the raw data that was outputted by any of the Simulated optics experiment/Simulator implementations. The ATS, Object Icon, Python, and RATFOR programs below can serve as reference implementations. The ATS program stores the raw data in a big array. The RATFOR does not store the raw data at all, but instead calculates as it inputs. The other programs use set data structures.
The analysis proceeds as follows.
- Read in the raw data. You do not have to retain the order of the data lines. In other words, you can, if you wish, store the data in a set structure.
- For each of the data lines, if its first entry is a "1", then swap the last two entries with each other, and do the same with the two entries before them. This is to account for the difference in geometry: the original light pulses are rotated 90° with respect to the apparatus, relative to if the first entry were a "0".
- Split the data into four groups, based on the values of the second and third line entries. In other words, one group for lines that look like , one for lines that look like , one for lines that look like , and one for lines that look like . Thus each group corresponds to a particular setting of the two polarizing beam splitters.
- Separately compute a correlation coefficient for each group. This is done as follows.
- Let equal the total number of lines in a group.
- Let equal the number of lines that have a "1" in the fourth entry. This is the number of detections by one of the light detectors on the left side.
- Let equal the number of lines that have a "1" in the fifth entry. This is the number of detections by the other light detector on the left side.
- Add up similar numbers for the right side: for the sixth entry in a line, and for the seventh entry.
- Compute the following numbers: , , , and . These numbers represent an estimate of the orientations of the two beam splitters, by inference from detection events in the light detectors.
- Compute and . These are angle difference formulas from trigonometry. We are "eliminating the coordinate system".
- Compute . This is the correlation coefficient for the group, estimated by inference from detection events.
- Designate each of the four correlation coefficients according to which group it is for. The Python example, for instance, refers to them as
kappaL1R1,kappaL1R2,kappaL2R1, andkappaL2R2. - Compute the estimated CHSH contrast:
-kappaL1R1 + kappaL1R2 + kappaL2R1 + kappaL2R2. - Display the number of "pulse pair events" that were simulated, the four correlation coefficients, and the estimated CHSH contrast. Also display the theoretical value of the CHSH contrast () and the difference between the estimated CHSH contrast and the theoretical value. Also display what we will call the CHSH violation: the difference between the estimated CHSH contrast and the number 2. Any value significantly greater than zero seems to be supposed to mean that the simulation is inherently "quantum mechanical" and involves "instantaneous action at a distance". But you can decide for yourself whether it really means that. Just do the calculations and print them. One should get an estimated CHSH contrast that is close to the theoretical value, and thus a relatively "huge" CHSH violation of about +0.8.
Extra "credit"
Here is a snapshot of what contains a thought experiment I have been working on: http://web.archive.org/web/20230529223105/https://crudfactory.com/Bell-assumes-his-conclusion.pdf The thought experiment concerns the same subject as this task, but more directly addresses the conundrum. It aims to illustrate the fallaciousness (despite their widespread acceptance by people regarded as "experts") of certain arguments underlying the conundrum, and thus to resolve the technical aspects of the conundrum. Whether it succeeds is up to you to decide. I am just a computer programmer, after all. And, even if I succeed in resolving the technical aspects of the conundrum, I have barely touched on the social aspects of it. Nevertheless, the thought experiment begs to be computer-animated, and I am not that kind of a programmer.
For extra "credit" (in the form of contribution to the community), write an animation or animations of the thought experiment.
--Chemoelectric (talk) 22:50, 29 May 2023 (UTC)
See also
- Simulated optics experiment/Simulator
- E. T. Jaynes, "Clearing up mysteries -- the original goal" Opening lecture of the eighth International MAXENT Workshop. 1988. This is a classic from the same physicist as the Jaynes-Cummings model is named for.
ATS
This implementation stores the data in an n-by-7 runtime-allocated array, and does not use set data structures. It should prove translatable to C, Ada, modern Fortran, etc., without much trouble.
I know ATS can be very difficult to comprehend if you do not know the language, but there is not much "ATSism" in this code. It likely is readable by many. The following peculiarities do occur: some of the loops are written as tail recursions (which Schemers, ML-users, etc., probably will have no difficulty with), and there is C code embedded in the program in more than one way. I embed C code in part to avoid using libraries other than the ATS prelude, which is minimal.
The program will have to be linked with the C math library and with an allocator, which can be malloc(3) (despite free(3) never being called), because the only allocation is for the array:
patscc -O2 -std=gnu2x -DATS_MEMALLOC_LIBC optics_simulation_analysis_task.dats -lm
(It is my habit to tell patscc to use C beyond -std=c99)
// Reference:
//
// A. F. Kracklauer, ‘EPR-B correlations: non-locality or geometry?’,
// J. Nonlinear Math. Phys. 11 (Supp.) 104–109 (2004).
// https://doi.org/10.2991/jnmp.2004.11.s1.13 (Open access, CC BY-NC)
#include "share/atspre_staload.hats"
typedef zero_or_one = intBtwe (0, 1)
val angleL1 = 0.0
and angleL2 = 45.0
and angleR1 = 22.5
and angleR2 = 67.5
(*------------------------------------------------------------------*)
%{^
#include <math.h>
#include <stdio.h>
//
// fscanf(3) will be a most convenient way to input the necessary
// data. It will ignore line breaks, but the line breaks are in the
// raw data for the convenience of line-by-line input, not because
// they matter.
//
static atstype_uint
_read_uint (FILE *f)
{
unsigned int u = 0U - 1U; // Yes, this is legal in standard C.
int bitbucket = fscanf (f, "%u", &u);
return u;
}
%}
fn
read_uint (inpf : FILEref) : uint =
let
typedef FILEstar = $extype"FILE *"
extern castfn FILEref2star : FILEref -<> FILEstar
extern fn _read_uint : FILEstar -> uint = "mac#_read_uint"
in
_read_uint (FILEref2star inpf)
end
fn
read_intGte0 (inpf : FILEref) : intGte 0 =
let
val u = g1ofg0 (read_uint inpf)
prval () = lemma_g1uint_param u (* Proves u >= 0. *)
in
u2i u
end
fn
read_zero_or_one (inpf : FILEref) : zero_or_one =
let
val i = read_intGte0 inpf
in
if i <= 1 then
i
else
begin
fprintln! (stderr_ref, "bad data");
exit (1)
end
end
(*------------------------------------------------------------------*)
fn
read_raw_data (inpf : FILEref) : mtrxszref zero_or_one =
let
val n = read_intGte0 inpf
val data = mtrxszref_make_elt<zero_or_one> (i2sz n, i2sz 7, 0)
var i : intGte 0
in
for (i := 0; i2sz i <> data.nrow(); i := succ i)
let
var j : intGte 0
in
for (j := 0; j <> 7; j := succ j)
data[i, j] := read_zero_or_one inpf
end;
data
end
fn {a : t@ype}
mtrxszref_swap (M : mtrxszref a,
i1 : int,
j1 : int,
i2 : int,
j2 : int) : void =
let
val m1 = M[i1, j1] and m2 = M[i2, j2]
in
M[i1, j1] := m2; M[i2, j2] := m1
end
fn
adjust_for_geometry (data : mtrxszref zero_or_one) : void =
let
var i : intGte 0
in
for (i := 0; i2sz i <> data.nrow(); i := succ i)
if data[i, 0] = 1 then
begin
mtrxszref_swap (data, i, 3, i, 4);
mtrxszref_swap (data, i, 5, i, 6)
end
end
fn
group_size (data : mtrxszref zero_or_one,
group : @(zero_or_one, zero_or_one)) : intGte 0 =
let
fun
loop (i : intGte 0,
count : intGte 0) : intGte 0 =
if i2sz i = data.nrow() then
count
else if data[i, 1] = group.0 && data[i, 2] = group.1 then
loop (succ i, succ count)
else
loop (succ i, count)
in
loop (0, 0)
end
fn
count_detections (data : mtrxszref zero_or_one,
group : @(zero_or_one, zero_or_one),
column : intBtwe (3, 6)) : intGte 0 =
let
fun
loop (i : intGte 0,
count : intGte 0) : intGte 0 =
if i2sz i = data.nrow() then
count
else if data[i, 1] = group.0 && data[i, 2] = group.1 then
loop (succ i, count + data[i, column])
else
loop (succ i, count)
in
loop (0, 0)
end
fn
correlation_coefficient (data : mtrxszref zero_or_one,
group : @(zero_or_one, zero_or_one))
: double =
let
extern fn sqrt : double -<> double = "mac#sqrt"
val N = group_size (data, group)
and NL1 = count_detections (data, group, 3)
and NL2 = count_detections (data, group, 4)
and NR1 = count_detections (data, group, 5)
and NR2 = count_detections (data, group, 6)
(* Equation (2.4) in the reference. *)
val sinL = sqrt (g0i2f NL1 / g0i2f N)
and cosL = sqrt (g0i2f NL2 / g0i2f N)
and sinR = sqrt (g0i2f NR1 / g0i2f N)
and cosR = sqrt (g0i2f NR2 / g0i2f N)
(* Equation (2.3). *)
val cosLR = (cosR * cosL) + (sinR * sinL)
and sinLR = (sinR * cosL) - (cosR * sinL)
(* Equation (2.5). *)
val kappa = (cosLR * cosLR) - (sinLR * sinLR)
in
kappa
end
fn
analyze_data (data : mtrxszref zero_or_one) : void =
let
extern fn sqrt : double -<> double = "mac#sqrt"
val kappa00 = correlation_coefficient (data, @(0, 0))
and kappa01 = correlation_coefficient (data, @(0, 1))
and kappa10 = correlation_coefficient (data, @(1, 0))
and kappa11 = correlation_coefficient (data, @(1, 1))
val chsh_contrast = ~kappa00 + kappa01 + kappa10 + kappa11
and chsh_contrast_nominal = 2.0 * sqrt 2.0
macdef X (x) = ignoret (,(x)) (* "ignore return value" *)
in
X ($extfcall (int, "printf", "\n light pulse events %9zu\n",
data.nrow()));
X ($extfcall (int, "printf", "\n correlation coefs\n"));
X ($extfcall (int, "printf", " %4.1f° %4.1f° %+9.6f\n",
angleL1, angleR1, kappa00));
X ($extfcall (int, "printf", " %4.1f° %4.1f° %+9.6f\n",
angleL1, angleR2, kappa01));
X ($extfcall (int, "printf", " %4.1f° %4.1f° %+9.6f\n",
angleL2, angleR1, kappa10));
X ($extfcall (int, "printf", " %4.1f° %4.1f° %+9.6f\n",
angleL2, angleR2, kappa11));
X ($extfcall (int, "printf", "\n CHSH contrast %+9.6f\n",
chsh_contrast));
X ($extfcall (int, "printf", " 2*sqrt(2) = nominal %+9.6f\n",
chsh_contrast_nominal));
X ($extfcall (int, "printf", " difference %+9.6f\n",
chsh_contrast - chsh_contrast_nominal));
X ($extfcall (int, "printf",
"\n CHSH violation %+9.6f\n\n",
chsh_contrast - 2.0))
end
implement main0 () =
let
val data = read_raw_data stdin_ref
in
adjust_for_geometry data;
analyze_data data
end- Output:
An example run:
light pulse events 100000
correlation coefs
0.0° 22.5° -0.707806
0.0° 67.5° +0.705415
45.0° 22.5° +0.712377
45.0° 67.5° +0.718882
CHSH contrast +2.844480
2*sqrt(2) = nominal +2.828427
difference +0.016053
CHSH violation +0.844480
Julia
Note that, statistically, the extra correlation in the analysis is inserted via the call to swap_LR_channels(), which is statistically dependent on whether the initial photon pair is [0, 90] or [90, 0].
using Printf
""" The data for one light pulse. """
mutable struct PulseData
logS::Int
logL::Int
logR::Int
detectedL1::Int
detectedL2::Int
detectedR1::Int
detectedR2::Int
end
""" Swap detectedL1 and detectedL2. """
swap_L_channels(p) = begin p.detectedL1, p.detectedL2 = p.detectedL2, p.detectedL1; p end
""" Swap detectedR1 and detectedR2. """
swap_R_channels(p) = begin p.detectedR1, p.detectedR2 = p.detectedR2, p.detectedR1; p end
"""
Swap channels on both right and left. This is done if the light source was (1,90°)
instead of (1,0°). For in that case the orientations of the polarizing beam splitters,
relative to the light source, is different by 90°.
"""
swap_LR_channels(p::PulseData) = (swap_L_channels(p); swap_R_channels(p))
"""
Split the data into two subsets, according to whether a set item satisfies a predicate.
The return value is a tuple, with those items satisfying the predicate in the first
tuple entry, the other items in the second entry.
"""
split_data(predicate, data) = (filter(predicate, data), filter(!predicate, data))
"""
Some data items are for a (1,0°) light pulse. The others are for a (1,90°) light pulse.
Thus the light pulses are oriented differently with respect to the polarizing beam splitters.
We adjust for that distinction here.
"""
function adjust_data_for_light_pulse_orientation(data)
data0, data90 = split_data(x -> x.logS == 0, data)
return vcat(data0, swap_LR_channels.(data90))
end
"""
Split the data into subsets per each arrangement of polarizing beam splitters.
"""
function split_data_according_to_PBS_setting(data)
dataL1, dataL2 = split_data(x -> x.logL == 0, data)
dataL1R1, dataL1R2 = split_data(x -> x.logR == 0, dataL1)
dataL2R1, dataL2R2 = split_data(x -> x.logR == 0, dataL2)
return (dataL1R1, dataL1R2, dataL2R1, dataL2R2)
end
"""
Compute the correlation coefficient for the subset of the data that corresponding
to a particular setting of the polarizing beam splitters.
"""
function compute_correlation_coefficient(angleL, angleR, data)
# We make certain the orientations of beam splitters are
# represented by perpendicular angles in the first and fourth
# quadrant. This restriction causes no loss of generality, because
# the orientation of the beam splitter is actually a rotated "X".
@assert all(0 <= x < 90 for x in (angleL, angleR)) "All orientations are acute angles"
# Note that the sine is non-negative in Quadrant 1, and the cosine
# is non-negative in Quadrant 4. Thus we can use the following
# estimates for cosine and sine. This is Equation (2.4) in the
# reference. (Note, one can also use Quadrants 1 and 2 and reverse
# the roles of cosine and sine. And so on like that.)
N = length(data)
NL1 = 0
NL2 = 0
NR1 = 0
NR2 = 0
for item in data
NL1 += item.detectedL1
NL2 += item.detectedL2
NR1 += item.detectedR1
NR2 += item.detectedR2
end
sinL = sqrt(NL1 / N)
cosL = sqrt(NL2 / N)
sinR = sqrt(NR1 / N)
cosR = sqrt(NR2 / N)
# Now we can apply the reference's Equation (2.3).
cosLR = (cosR * cosL) + (sinR * sinL)
sinLR = (sinR * cosL) - (cosR * sinL)
# And then Equation (2.5).
kappa = (cosLR * cosLR) - (sinLR * sinLR)
return kappa
end
""" Read the raw data. Its order does not actually matter. """
function read_raw_data(stream)
return [PulseData(row...) for row in filter(a -> length(a) == 7,
[[parse(Int, s) for s in split(line)] for line in readlines(stream)])]
end
function run_analysis(inputstream)
# Polarizing beam splitter orientations commonly used in actual
# experiments. These must match the values used in the simulation
# itself. They are by design all either zero degrees or in the
# first quadrant.
anglesL = [0.0, 45.0]
anglesR = [22.5, 67.5]
@assert all(0 <= x < 90 for x in [anglesL; anglesR])
data = read_raw_data(inputstream) |> adjust_data_for_light_pulse_orientation
dataL1R1, dataL1R2, dataL2R1, dataL2R2 = split_data_according_to_PBS_setting(data)
kappaL1R1 = compute_correlation_coefficient(anglesL[1], anglesR[1], dataL1R1)
kappaL1R2 = compute_correlation_coefficient(anglesL[1], anglesR[2], dataL1R2)
kappaL2R1 = compute_correlation_coefficient(anglesL[2], anglesR[1], dataL2R1)
kappaL2R2 = compute_correlation_coefficient(anglesL[2], anglesR[2], dataL2R2)
chsh_contrast = -kappaL1R1 + kappaL1R2 + kappaL2R1 + kappaL2R2
# The nominal value of the CHSH contrast for the chosen polarizer
# orientations is 2*sqrt(2).
chsh_contrast_nominal = 2 * sqrt(2.0)
@printf("\n light pulse events %9d\n\n", length(data))
println(" correlation coefs")
@printf(" %2d° %2d° %+9.6f\n", anglesL[1], anglesR[1], kappaL1R1)
@printf(" %2d° %2d° %+9.6f\n", anglesL[1], anglesR[2], kappaL1R2)
@printf(" %2d° %2d° %+9.6f\n", anglesL[2], anglesR[1], kappaL2R1)
@printf(" %2d° %2d° %+9.6f\n\n", anglesL[2], anglesR[2], kappaL2R2)
@printf(" CHSH contrast %+9.6f\n", chsh_contrast)
@printf(" 2*sqrt(2) = nominal %+9.6f\n", chsh_contrast_nominal)
@printf(" difference %+9.6f\n", chsh_contrast - chsh_contrast_nominal)
# A "CHSH violation" occurs if the CHSH contrast is >2.
# https://en.wikipedia.org/w/index.php?title=CHSH_inequality&oldid=1142431418
@printf("\n CHSH violation %+9.6f\n\n", chsh_contrast - 2)
end
run_analysis("datalog.log")
- Output:
An example using data from the Julia implementation of Simulated optics experiment/Simulator.
light pulse events 1000000
correlation coefs
0° 22° -0.459472
0° 68° +0.955042
45° 22° +0.752309
45° 68° +0.749286
CHSH contrast +2.913088
2*sqrt(2) = nominal +2.828427
difference +0.084660
CHSH violation +0.913088
Nim
import std/[hashes, math, os, sequtils, sets, strscans, strformat, strutils, sugar]
type
PulseData = ref object
logS: int
logL: int
logR: int
detectedL1: int
detectedL2: int
detectedR1: int
detectedR2: int
PulseDataSet = HashSet[PulseData]
proc hash(pd: PulseData): Hash =
## Hash value of a PulseData object (needed for HashSet).
hash(cast[int](pd[].addr))
func swapLRChannels(pd: PulseData) =
## Swap channels on both right and left. This is done if the
## light source was (1,90°) instead of (1,0°). For in that case
## the orientations of the polarizing beam splitters, relative to
## the light source, is different by 90°.
swap pd.detectedL1, pd.detectedL2
swap pd.detectedR1, pd.detectedR2
proc split(data: PulseDataSet; predicate: proc(pd: PulseData): bool): tuple[a, b: PulseDataSet] =
## Split the data into two subsets, according to whether a set
## item satisfies a predicate. The return value is a tuple, with
## those items satisfying the predicate in the first tuple entry, the
## other items in the second entry.
let subset1 = collect:
for x in data:
if predicate(x): {x}
let subset2 = data - subset1
result = (subset1, subset2)
proc adjustForLightPulseOrientation(data: PulseDataSet): PulseDataSet =
## Some data items are for a (1,0°) light pulse. The others are
## for a (1,90°) light pulse. Thus the light pulses are oriented
## differently with respect to the polarizing beam splitters. We
## adjust for that distinction here.
let (data0, data90) = data.split(proc(pd: PulseData): bool = pd.logS == 0)
for item in data90:
item.swapLRChannels()
result = union(data0, data90)
proc splitAccordingToPbsSetting(data: PulseDataSet): tuple[a, b, c, d: PulseDataSet] =
## Split the data into four subsets: one subset for each
## arrangement of the two polarizing beam splitters.
let (dataL1, dataL2) = data.split(proc(pd: PulseData): bool = pd.logL == 0)
let (dataL1R1, dataL1R2) = dataL1.split(proc(pd: PulseData): bool = pd.logR == 0)
let (dataL2R1, dataL2R2) = dataL2.split(proc(pd: PulseData): bool = pd.logR == 0)
result = (dataL1R1, dataL1R2, dataL2R1, dataL2R2)
func computeCorrelationCoefficient(data: PulseDataSet; angleL, angleR: float): float =
## Compute the correlation coefficient for the subset of the data that
## corresponding to a particular setting of the polarizing beam splitters.
# We make certain the orientations of beam splitters are
# represented by perpendicular angles in the first and fourth
# quadrant. This restriction causes no loss of generality, because
# the orientation of the beam splitter is actually a rotated "X".
assert angleL >= 0 and angleL < 90 and angleR >= 0 and angleR < 90
#perpendicularL = angleL - 90 # In Quadrant 4.
#perpendicularR = angleR - 90 # In Quadrant 4.
# Note that the sine is non-negative in Quadrant 1, and the cosine
# is non-negative in Quadrant 4. Thus we can use the following
# estimates for cosine and sine. This is Equation (2.4) in the
# reference. (Note, one can also use Quadrants 1 and 2 and reverse
# the roles of cosine and sine. And so on like that.)
let n = data.len
var nL1, nL2, nR1, nR2 = 0
for item in data:
inc nL1, item.detectedL1
inc nL2, item.detectedL2
inc nR1, item.detectedR1
inc nR2, item.detectedR2
let sinL = sqrt(nL1 / n)
let cosL = sqrt(nL2 / n)
let sinR = sqrt(nR1 / n)
let cosR = sqrt(nR2 / n)
# Now we can apply the reference's Equation (2.3).
let cosLR = (cosR * cosL) + (sinR * sinL)
let sinLR = (sinR * cosL) - (cosR * sinL)
# And then Equation (2.5).
result = (cosLR * cosLR) - (sinLR * sinLR)
proc readRawData(filename: string): PulseDataSet =
## Read the raw data. Its order does not actually matter, so we
## return the data as a HashSet.
func makeRecord(line: string): PulseData =
new(result)
discard line.scanf("$i $i $i $i $i $i $i", result.logS, result.logL, result.logR,
result.detectedL1, result.detectedL2, result.detectedR1, result.detectedR2)
proc readData(f: File): PulseDataSet =
let numPulses = f.readline().parseInt()
for i in 1..numPulses:
result.incl f.readLine().makeRecord()
if filename != "-":
let f = open(filename, fmRead)
result = f.readData()
f.close()
else:
result = stdin.readData()
when isMainModule:
if paramCount() notin [0, 1]:
quit &"Usage: {getAppFilename()} [RAW_DATA_FILE]", QuitFailure
let rawDataFilename = if paramCount() == 1: paramStr(1) else: "-"
# Polarizing beam splitter orientations commonly used in actual
# experiments. These must match the values used in the simulation
# itself. They are by design all either zero degrees or in the
# first quadrant.
const AnglesL = [0.0, 45.0]
const AnglesR = [22.5, 67.5]
assert allIt(AnglesL, it >= 0 and it < 90) and allIt(AnglesR, it >= 0 and it < 90)
var data = readRawData(rawDataFilename)
data = adjustForLightPulseOrientation(data)
let (dataL1R1, dataL1R2, dataL2R1, dataL2R2) = data.splitAccordingToPbsSetting()
let
kappaL1R1 = dataL1R1.computeCorrelationCoefficient(AnglesL[0], AnglesR[0])
kappaL1R2 = dataL1R2.computeCorrelationCoefficient(AnglesL[0], AnglesR[1])
kappaL2R1 = dataL2R1.computeCorrelationCoefficient(AnglesL[1], AnglesR[0])
kappaL2R2 = dataL2R2.computeCorrelationCoefficient(AnglesL[1], AnglesR[1])
let chshContrast = -kappaL1R1 + kappaL1R2 + kappaL2R1 + kappaL2R2
# The nominal value of the CHSH contrast for the chosen polarizer
# orientations is 2*sqrt(2).
let chshContrastNominal = 2 * sqrt(2.0)
echo()
echo &" light pulse events {data.len:9}"
echo()
echo " correlation coefs"
echo &" {AnglesL[0]:4.1f}° {AnglesR[0]:4.1f}° {kappaL1R1:+9.6f}"
echo &" {AnglesL[0]:4.1f}° {AnglesR[1]:4.1f}° {kappaL1R2:+9.6f}"
echo &" {AnglesL[1]:4.1f}° {AnglesR[0]:4.1f}° {kappaL2R1:+9.6f}"
echo &" {AnglesL[1]:4.1f}° {AnglesR[1]:4.1f}° {kappaL2R2:+9.6f}"
echo()
echo &" CHSH contrast {chshContrast:+9.6f}"
echo &" 2*sqrt(2) = nominal {chshContrastNominal:+9.6f}"
echo &" difference {chshContrast - chshContrastNominal:+9.6f}"
# A "CHSH violation" occurs if the CHSH contrast is >2.
# https://en.wikipedia.org/w/index.php?title=CHSH_inequality&oldid=1142431418
echo()
echo &" CHSH violation {chshContrast - 2:+9.6f}"
echo()
- Output:
The result is identical to that of the Python version. For instance, using data produced by the Python version of the simulator:
light pulse events 100000
correlation coefs
0.0° 22.5° -0.708395
0.0° 67.5° +0.705963
45.0° 22.5° +0.709600
45.0° 67.5° +0.705349
CHSH contrast +2.829306
2*sqrt(2) = nominal +2.828427
difference +0.000879
CHSH violation +0.829306
ObjectIcon
#!/bin/env -S oiscript
#
# If you run this program as a script, you might get something like
# the following as a "help" message:
#
# Usage: /tmp/oixcache/f8f91465ab42189d5623040505fa7f12 [-|raw_data_file]
#
# But you can instead use "oit" to compile this program into a
# bytecode executable, which would print this sort of message instead:
#
# Usage: ./my_program [-|raw_data_file]
#
#
# Reference:
#
# A. F. Kracklauer, ‘EPR-B correlations: non-locality or geometry?’,
# J. Nonlinear Math. Phys. 11 (Supp.) 104–109 (2004).
# https://doi.org/10.2991/jnmp.2004.11.s1.13 (Open access, CC BY-NC)
#
# This is necessary because there are "degree" symbols (°) in the
# code:
$encoding UTF-8
# To get the following more specific set of imports, one can write
# "import io, util" and then run "fiximports" on the source code:
import
io(BufferStream, FileStream, open, stop),
ipl.printf(printf),
util(Math)
$define angleL1 0.0
$define angleL2 45.0
$define angleR1 22.5
$define angleR2 67.5
procedure main (args)
# "&why" does not exist in regular Icon, although "stop"
# does.
#
# Fellow Gentoo users might notice an uncanny resemblance between
# the use of "stop" here and the use of "die" in Gentoo ebuilds. In
# both cases, if we left out the check, the error might be silently
# ignored.
case *args of
{
0 : sally_forth ("-") | stop (&why)
1 : sally_forth (args[1]) | stop (&why)
default :
{
io.write ("Usage: ", &progname, " [-|raw_data_file]")
exit (1)
}
}
end
procedure sally_forth (input_file)
local data, dataL1, dataL2, dataR1, dataR2
local dataL1R1, dataL1R2, dataL2R1, dataL2R2
local kappaL1R1, kappaL1R2, kappaL2R1, kappaL2R2
local chsh_contrast, chsh_contrast_nominal
data := read_raw_data (input_file)
data := adjust_for_geometry (data)
# The following means of dividing the set into four non-overlapping
# groups has been designed for clarity rather than efficiency. (We
# should not care about efficiency, anyway, in a program such as
# this.)
dataL1 := filter_set (logL_equals_0, data)
dataL2 := filter_set (logL_equals_1, data)
dataR1 := filter_set (logR_equals_0, data)
dataR2 := filter_set (logR_equals_1, data)
dataL1R1 := dataL1 ** dataR1 # Set intersection.
dataL1R2 := dataL1 ** dataR2
dataL2R1 := dataL2 ** dataR1
dataL2R2 := dataL2 ** dataR2
# The sizes of the groups should add up to the size of the entire
# data.
*dataL1R1 + *dataL1R2 + *dataL2R1 + *dataL2R2 = *data |
{
#
# "&why" does not exist in regular Icon.
#
# In all the dialects of Icon, however, it is common for failure
# of a procedure to indicate an error.
#
# This has a severe downside that, in many contexts, a failure
# is silently ignored. This behavior is unlike an exception
# (which is either caught or it terminates the
# program). However, in some ways it is like the practice of
# indicating an error through a return value. Differences from
# that include: (a) failure will often cause the containing
# context to fail as well (possibly resulting in assignments and
# swaps/exchanges being undone!), and (b) there is no return
# value in case of failure.
#
# BTW Object Icon does have exceptions:
# https://chemoelectric.github.io/objecticon/ExceptionPackage.html
#
&why := "The four groups of data have been formed incorrectly."
fail
}
kappaL1R1 := correlation_coefficient (angleL1, angleR1, dataL1R1)
kappaL1R2 := correlation_coefficient (angleL1, angleR2, dataL1R2)
kappaL2R1 := correlation_coefficient (angleL2, angleR1, dataL2R1)
kappaL2R2 := correlation_coefficient (angleL2, angleR2, dataL2R2)
chsh_contrast := (-kappaL1R1) + kappaL1R2 + kappaL2R1 + kappaL2R2
printf ("\n light pulse events %9d\n", *data)
printf ("\n correlation coefs\n")
# The u"" is necessary in Object Icon when the string is a Unicode
# string.
printf (u" %4.1r° %4.1r° %+9.6r\n",
angleL1, angleR1, kappaL1R1)
printf (u" %4.1r° %4.1r° %+9.6r\n",
angleL1, angleR2, kappaL1R2)
printf (u" %4.1r° %4.1r° %+9.6r\n",
angleL2, angleR1, kappaL2R1)
printf (u" %4.1r° %4.1r° %+9.6r\n",
angleL2, angleR2, kappaL2R2)
chsh_contrast_nominal := 2.0 * Math.sqrt (2.0)
printf ("\n CHSH contrast %+9.6r\n", chsh_contrast)
printf (" 2*sqrt(2) = nominal %+9.6r\n", chsh_contrast_nominal)
printf (" difference %+9.6r\n", (chsh_contrast -
chsh_contrast_nominal))
# A "CHSH violation" occurs if the CHSH contrast is >2.
# https://en.wikipedia.org/w/index.php?title=CHSH_inequality&oldid=1142431418
#
# Note CAREFULLY, however, that we are talking here about simulation
# of a CLASSICAL optics model, using the CLASSICAL formula for
# correlation coefficients (of the respective electromagnetic
# fields). We are calculating those coefficients by inference from
# detection events.
#
# Some might try to object that a CHSH contrast cannot be calculated
# from correlation coefficients computed according to the classical
# formula. But that would be assuming the conclusion: OF COURSE no
# classical model can have a CHSH violation, if by fiat you do not
# let a classical model and the classical formulas be used!
#
printf ("\n CHSH violation %+9.6r\n\n", (chsh_contrast - 2))
# The following is necessary (it is equivalent to "return &null"),
# because we are using failure to indicate an error.
return
end
procedure correlation_coefficient (angleL, angleR, data_group)
local N, NL1, NL2, NR1, NR2, item
local cosL, sinL, cosR, sinR, cosLR, sinLR, kappa
# We know angleL and angleR are in Quadrant 1. Its two
# perpendiculars are thus in Quadrants 2 and 4. We can choose
# whichever of the two we want. Choose Quadrant 4. The sine is
# non-negative in Quadrant 1, and the cosine is non-negative in
# Quadrant 4. Thus we can use the following estimates for cosine and
# sine, and not have to deal with negative signs.
#
# This is Equation (2.4) in the reference.
#
N := real (*data_group)
NL1 := 0.0
NL2 := 0.0
NR1 := 0.0
NR2 := 0.0
every item := !data_group do
{
NL1 +:= item.detectedL1()
NL2 +:= item.detectedL2()
NR1 +:= item.detectedR1()
NR2 +:= item.detectedR2()
}
sinL := Math.sqrt (NL1 / N)
cosL := Math.sqrt (NL2 / N)
sinR := Math.sqrt (NR1 / N)
cosR := Math.sqrt (NR2 / N)
# Apply the reference's Equation (2.3).
cosLR := (cosR * cosL) + (sinR * sinL)
sinLR := (sinR * cosL) - (cosR * sinL)
# Apply Equation (2.5).
kappa := (cosLR * cosLR) - (sinLR * sinLR)
return kappa
end
procedure read_raw_data (input_file)
local data, f
local num_lines
local logS, logL, logR
local detectedL1, detectedL2
local detectedR1, detectedR2
if type (input_file) ~== "string" then
{
f := input_file
num_lines := integer (f.read())
data := set ()
# The "every" construct below goes through all the numbers
# generated by "1 to num_lines", after which it quits because the
# generator FAILS. (I leave out any loop variable, "every i :=",
# because one is not needed here.) One might, especially at first,
# accidentally write something like "while i := 1 to num_lines",
# but this would actually be a forever-loop, with i always equal
# to 1. The "while" construct, unlike "every", would restart the
# generator from scratch each time, so that the generator ALWAYS
# SUCCEEDS. But "while" is what you would want for something like
# "while line := f.read()", because "read" is not a generator. It
# is a procedure that either returns a line of input or fails.
every 1 to num_lines do
{
# The following is a "string-scanning environment". It is a
# series of operations on the string returned by f.read(). See
# also
# https://en.wikipedia.org/w/index.php?title=Icon_(programming_language)&oldid=1146573004
#
# Object Icon is not mentioned in that article, but right now I
# do not feel like editing the Wikipedia to mention it.
#
# It IS mentioned that Python generators are supposedly inspired
# by Icon generators (and presumably also by co-expressions). I
# do not know if that is true, but there are similarities:
# Python "yield" is similar to Icon "suspend". (Scheme’s
# "call/cc" is far more general, however, and can be used to
# implement yield/suspend in a few lines of code.)
f.read() ?
{
tab(many(' \t'))
logS := integer (tab(many(&digits)))
tab(many(' \t'))
logL := integer (tab(many(&digits)))
tab(many(' \t'))
logR := integer (tab(many(&digits)))
tab(many(' \t'))
detectedL1 := integer (tab(many(&digits)))
tab(many(' \t'))
detectedL2 := integer (tab(many(&digits)))
tab(many(' \t'))
detectedR1 := integer (tab(many(&digits)))
tab(many(' \t'))
detectedR2 := integer (tab(many(&digits)))
}
insert (data, PulseEventData (logS, logL, logR,
detectedL1, detectedL2,
detectedR1, detectedR2))
}
}
else if input_file == "-" then
data := read_raw_data (BufferStream (FileStream.stdin))
else
{
f := open (input_file, "r") | stop (&why)
data := read_raw_data (f) | stop (&why)
f.close ()
}
return data
end
procedure adjust_for_geometry (data)
local S, S1, S2, S2a, item
S := data
S1 := filter_set (logS_equals_0, S)
S2 := S -- S1
S2a := set ()
every item := !S2 do
#
# Swap channels, to account for geometry.
#
# This step has been objected to as introducing "quantum magic",
# but note that we are doing a CLASSICAL data analysis, having
# nought to do with quantum mechanics. The light pulse source has
# recorded for us which type of light pulses were emitted. One
# might or might not be able to imagine a device for doing this
# with "photons", but we do not need to: (a) this is a SIMULATION,
# so we CAN use information we do not KNOW HOW to get, as long as
# the information is moving strictly forwards in time; besides
# which, (b) we can INSTEAD imagine a device that randomly selects
# from one of two pairs of light beams, and which forms pulses by
# opening and closing a shutter.
#
# According to "Bell's Theorem", it is impossible for such a
# classical model--one that works solely by contact action moving
# forwards in time--to produce the results we get. In my opinion,
# this simulation should be enough to end the entire debacle this
# "theorem" has helped create. But, in fact, it is easy to show
# that John Bell's "proof" does the equivalent of using the
# Angle-Side-Side "theorem" in geometry class. It is fake
# mathematics. (See the "Talk" page of the task, where I
# explain Bell's mistake in just the few words it takes.)
#
# One might or might not have noticed we are also including
# non-detection events in denominators, even though they are not
# part of the measured data set, in experiments that do not record
# when pulses were emitted in the first place. This is a
# SIMULATION, however, so we CAN include them.
#
# -- Chemoelectric, 30 May 2023
#
# In any case, we are demonstrating how to use sets in Object
# Icon. Sets are very similar in regular Icon and in Unicon,
# although regular Icon does not have "classes". (It DOES have
# "records".)
#
insert (S2a, PulseEventData (item.logS(),
item.logL(), item.logR(),
item.detectedL2(),
item.detectedL1(),
item.detectedR2(),
item.detectedR1()))
return S1 ++ S2a
end
procedure filter_set (semidet_predicate, S)
# The predicate is "semi-deterministic". It does not return a
# boolean. Instead, it either succeeds (and in Object Icon has to
# return something, which we ignore) or it fails. You will see the
# same approach to decision-making in Prolog and Mercury programs.
local S1, item
S1 := set ()
every item := !S do
if semidet_predicate (item) then
insert (S1, item)
return S1
end
procedure logS_equals_0 (event_data)
# This is a semidet predicate. It succeeds if logS is zero,
# otherwise it fails. The return value (which is returned only on
# success and is always 0) will be ignored.
return event_data.logS() = 0
end
#######################################
# Some more semidet predicates:
procedure logL_equals_0 (event_data)
return event_data.logL() = 0
end
procedure logL_equals_1 (event_data)
return event_data.logL() = 1
end
procedure logR_equals_0 (event_data)
return event_data.logR() = 0
end
procedure logR_equals_1 (event_data)
return event_data.logR() = 1
end
#######################################
final class PulseEventData ()
# I have tried here to make PulseEventData objects both immutable
# and unextendable, so no one can doubt the absence of "clever
# tricks". They are merely containers of data: each of them unique,
# so they can be stored in sets.
private var_logS, var_logL, var_logR
private var_detectedL1, var_detectedL2
private var_detectedR1, var_detectedR2
public new (logS, logL, logR,
detectedL1, detectedL2,
detectedR1, detectedR2)
var_logS := logS
var_logL := logL
var_logR := logR
var_detectedL1 := detectedL1
var_detectedL2 := detectedL2
var_detectedR1 := detectedR1
var_detectedR2 := detectedR2
return
end
public logS ()
return var_logS
end
public logL ()
return var_logL
end
public logR ()
return var_logR
end
public detectedL1 ()
return var_detectedL1
end
public detectedL2 ()
return var_detectedL2
end
public detectedR1 ()
return var_detectedR1
end
public detectedR2 ()
return var_detectedR2
end
public to_string ()
return ("PulseEventData(" || logS() ||
", " || logL() ||
", " || logR() ||
", " || detectedL1() ||
", " || detectedL2() ||
", " || detectedR1() ||
", " || detectedR2() || ")")
end
end- Output:
Using data generated by a 10000-point run of the Python version of the simulation.
Please be aware that one run of that Python program will not produce the same output as another. There are two reasons: (1) the random number generator is seeded randomly, and (2) the simulated devices are running as separate processes. (How Python handles random numbers in a program running as multiple processes I haven't the vaguest idea.)
A standard set of raw data might be published for people to use. Then we could always get the same numbers. But, to quote Cyrano Jones, "What would happen to Man's quest for knowledge?"
light pulse events 10000
correlation coefs
0.0° 22.5° -0.699476
0.0° 67.5° +0.714172
45.0° 22.5° +0.692368
45.0° 67.5° +0.711478
CHSH contrast +2.817493
2*sqrt(2) = nominal +2.828427
difference -0.010935
CHSH violation +0.817493
Phix
Gone for as simple as possible, includes the generation phase, and for javascript compatability file save/load commented out.
with javascript_semantics sequence datalog enum S, L, R, L1, L2, R1, R2, LL=$ constant angles = {0, -- (this is indexed via L and R) {0.0,45.0}, {22.5,67.5}} procedure detector_recieve(integer lr, i, atom amplitude) atom intensity := amplitude * amplitude, randnum := rnd() datalog[i][lr] := (randnum <= intensity) end procedure procedure splitter_receive(integer lr, i, atom source_angle) integer rand_0_or_1 := rand_range(0,1) datalog[i][lr] = rand_0_or_1 atom my_angle := angles[lr][rand_0_or_1+1], relative_angle := my_angle - source_angle, amplitude1 := abs(cos(PI/180 * relative_angle)), amplitude2 := abs(sin(PI/180 * relative_angle)) detector_recieve(2*lr,i,amplitude1) detector_recieve(2*lr+1,i,amplitude2) end procedure procedure light_source(integer num_events) for i=1 to num_events do integer rand_0_or_1 := rand_range(0,1) datalog[i][S] = rand_0_or_1 if rand_0_or_1 = 0 then splitter_receive(L,i,0.0) splitter_receive(R,i,90.0) else splitter_receive(L,i,90.0) splitter_receive(R,i,0.0) end if end for end procedure procedure generate(integer num_events, object output_file=1) datalog = repeat(repeat(-9,LL),num_events) light_source(num_events) if string(output_file) or num_events<=20 then integer fn = iff(string(output_file)?open(output_file,"w"):output_file) printf(fn,"%d\n",num_events) for l in datalog do printf(fn,"%s\n",join(l,fmt:="%d")) end for if string(output_file) then close(fn) end if end if end procedure generate(100000) --generate(1e5,"datalog.log") function adjust_data_for_light_pulse_orientation(sequence data) for i,d in data do if d[S]=1 then data[i] = reinstate(d,{L1,L2,R1,R2},extract(d,{L2,L1,R2,R1})) end if end for return data end function function split_data_according_to_PBS_setting(sequence data) sequence dataL1R1 = {}, dataL1R2 = {}, dataL2R1 = {}, dataL2R2 = {} for d in data do sequence dLR = d[L..R] if dLR={0,0} then dataL1R1 = append(dataL1R1,d) elsif dLR={0,1} then dataL1R2 = append(dataL1R2,d) elsif dLR={1,0} then dataL2R1 = append(dataL2R1,d) elsif dLR={1,1} then dataL2R2 = append(dataL2R2,d) else ?9/0 end if end for return {dataL1R1, dataL1R2, dataL2R1, dataL2R2} end function function compute_correlation_coefficient(sequence data) integer N = length(data), NL1 = 0, NL2 = 0, NR1 = 0, NR2 = 0 for item in data do NL1 += item[L1] NL2 += item[L2] NR1 += item[R1] NR2 += item[R2] end for atom sinL = sqrt(NL1 / N), cosL = sqrt(NL2 / N), sinR = sqrt(NR1 / N), cosR = sqrt(NR2 / N), cosLR = (cosR * cosL) + (sinR * sinL), sinLR = (sinR * cosL) - (cosR * sinL), kappa = (cosLR * cosLR) - (sinLR * sinLR) return kappa end function --function read_from_file(string input_file) -- sequence data = read_lines(input_file) -- assert(to_integer(data[1])+1=length(data)) -- data = data[2..$] -- for i,d in data do -- data[i] = scanf(d,"%d %d %d %d %d %d %d")[1] -- end for -- return data --end function procedure run_analysis() --procedure run_analysis(string input_file) -- sequence data = read_from_file(input_file) -- assert(data=datalog) sequence data = datalog datalog = {} -- (kill refcount) data = adjust_data_for_light_pulse_orientation(data) sequence {dataL1R1, dataL1R2, dataL2R1, dataL2R2} = split_data_according_to_PBS_setting(data) atom {{},{angleL1,angleL2},{angleR1,angleR2}} = angles, kappaL1R1 = compute_correlation_coefficient(dataL1R1), kappaL1R2 = compute_correlation_coefficient(dataL1R2), kappaL2R1 = compute_correlation_coefficient(dataL2R1), kappaL2R2 = compute_correlation_coefficient(dataL2R2), chsh_contrast = -kappaL1R1 + kappaL1R2 + kappaL2R1 + kappaL2R2, chsh_contrast_nominal = 2 * sqrt(2.0) printf(1,"\n light pulse events %9d\n\n", length(data)) printf(1," correlation coefs\n") printf(1," %4.1f° %4.1f° %+9.6f\n", {angleL1, angleR1, kappaL1R1}) printf(1," %4.1f° %4.1f° %+9.6f\n", {angleL1, angleR2, kappaL1R2}) printf(1," %4.1f° %4.1f° %+9.6f\n", {angleL2, angleR1, kappaL2R1}) printf(1," %4.1f° %4.1f° %+9.6f\n", {angleL2, angleR2, kappaL2R2}) printf(1,"\n") printf(1," CHSH contrast %+9.6f\n", chsh_contrast) printf(1," 2*sqrt(2) = nominal %+9.6f\n", chsh_contrast_nominal) printf(1," difference %+9.6f\n", chsh_contrast - chsh_contrast_nominal) printf(1,"\n CHSH violation %+9.6f\n\n", chsh_contrast - 2) end procedure run_analysis() --run_analysis("datalog.log")
- Output:
light pulse events 100000 correlation coefs 0.0° 22.5° -0.708021 0.0° 67.5° +0.700617 45.0° 22.5° +0.709869 45.0° 67.5° +0.712739 CHSH contrast +2.831245 2*sqrt(2) = nominal +2.828427 difference +0.002818 CHSH violation +0.831245
Python
The program takes one optional command line argument: the filename of the raw data saved by one of the Simulated optics experiment/Simulator implementations. If the argument is omitted or is "-", the raw data is read from standard input.
#!/bin/env -S python3
#
# Reference:
#
# A. F. Kracklauer, ‘EPR-B correlations: non-locality or geometry?’,
# J. Nonlinear Math. Phys. 11 (Supp.) 104–109 (2004).
# https://doi.org/10.2991/jnmp.2004.11.s1.13 (Open access, CC BY-NC)
#
import sys
from math import atan2, cos, floor, pi, radians, sin, sqrt
class PulseData():
'''The data for one light pulse.'''
def __init__(self, logS, logL, logR, detectedL1, detectedL2,
detectedR1, detectedR2):
self.logS = logS
self.logL = logL
self.logR = logR
self.detectedL1 = detectedL1
self.detectedL2 = detectedL2
self.detectedR1 = detectedR1
self.detectedR2 = detectedR2
def __repr__(self):
return ("PulseData(" + str(self.logS) + ", "
+ str(self.logL) + ", "
+ str(self.logR) + ", "
+ str(self.detectedL1) + ", "
+ str(self.detectedL2) + ", "
+ str(self.detectedR1) + ", "
+ str(self.detectedR2) + ")")
def swap_L_channels(self):
'''Swap detectedL1 and detectedL2.'''
(self.detectedL1, self.detectedL2) = \
(self.detectedL2, self.detectedL1)
def swap_R_channels(self):
'''Swap detectedR1 and detectedR2.'''
(self.detectedR1, self.detectedR2) = \
(self.detectedR2, self.detectedR1)
def swap_LR_channels(self):
'''Swap channels on both right and left. This is done if the
light source was (1,90°) instead of (1,0°). For in that case
the orientations of the polarizing beam splitters, relative to
the light source, is different by 90°.'''
self.swap_L_channels()
self.swap_R_channels()
def split_data(predicate, data):
'''Split the data into two subsets, according to whether a set
item satisfies a predicate. The return value is a tuple, with
those items satisfying the predicate in the first tuple entry, the
other items in the second entry.
'''
subset1 = {x for x in data if predicate(x)}
subset2 = data - subset1
return (subset1, subset2)
def adjust_data_for_light_pulse_orientation(data):
'''Some data items are for a (1,0°) light pulse. The others are
for a (1,90°) light pulse. Thus the light pulses are oriented
differently with respect to the polarizing beam splitters. We
adjust for that distinction here.'''
data0, data90 = split_data(lambda item: item.logS == 0, data)
for item in data90:
item.swap_LR_channels()
return (data0 | data90) # Set union.
def split_data_according_to_PBS_setting(data):
'''Split the data into four subsets: one subset for each
arrangement of the two polarizing beam splitters.'''
dataL1, dataL2 = split_data(lambda item: item.logL == 0, data)
dataL1R1, dataL1R2 = \
split_data(lambda item: item.logR == 0, dataL1)
dataL2R1, dataL2R2 = \
split_data(lambda item: item.logR == 0, dataL2)
return (dataL1R1, dataL1R2, dataL2R1, dataL2R2)
def compute_correlation_coefficient(angleL, angleR, data):
'''Compute the correlation coefficient for the subset of the data
that corresponding to a particular setting of the polarizing beam
splitters.'''
# We make certain the orientations of beam splitters are
# represented by perpendicular angles in the first and fourth
# quadrant. This restriction causes no loss of generality, because
# the orientation of the beam splitter is actually a rotated "X".
assert (all(0 <= x and x < 90 for x in (angleL, angleR)))
#perpendicularL = angleL - 90 # In Quadrant 4.
#perpendicularR = angleR - 90 # In Quadrant 4.
# Note that the sine is non-negative in Quadrant 1, and the cosine
# is non-negative in Quadrant 4. Thus we can use the following
# estimates for cosine and sine. This is Equation (2.4) in the
# reference. (Note, one can also use Quadrants 1 and 2 and reverse
# the roles of cosine and sine. And so on like that.)
N = len(data)
NL1 = 0
NL2 = 0
NR1 = 0
NR2 = 0
for item in data:
NL1 += item.detectedL1
NL2 += item.detectedL2
NR1 += item.detectedR1
NR2 += item.detectedR2
sinL = sqrt(NL1 / N)
cosL = sqrt(NL2 / N)
sinR = sqrt(NR1 / N)
cosR = sqrt(NR2 / N)
# Now we can apply the reference's Equation (2.3).
cosLR = (cosR * cosL) + (sinR * sinL)
sinLR = (sinR * cosL) - (cosR * sinL)
# And then Equation (2.5).
kappa = (cosLR * cosLR) - (sinLR * sinLR)
return kappa
def read_raw_data(filename):
'''Read the raw data. Its order does not actually matter, so we
return the data as a set.'''
def make_record(line):
x = line.split()
record = PulseData(logS = int(x[0]),
logL = int(x[1]),
logR = int(x[2]),
detectedL1 = int(x[3]),
detectedL2 = int(x[4]),
detectedR1 = int(x[5]),
detectedR2 = int(x[6]))
return record
def read_data(f):
num_pulses = int(f.readline())
data = set()
for i in range(num_pulses):
data.add(make_record(f.readline()))
return data
if filename != "-":
with open(filename, "r") as f:
data = read_data(f)
else:
data = read_data(sys.stdin)
return data
if __name__ == '__main__':
if len(sys.argv) != 1 and len(sys.argv) != 2:
print("Usage: " + sys.argv[0] + " [RAW_DATA_FILE]")
sys.exit(1)
raw_data_filename = (sys.argv[1] if len(sys.argv) == 2 else "-")
# Polarizing beam splitter orientations commonly used in actual
# experiments. These must match the values used in the simulation
# itself. They are by design all either zero degrees or in the
# first quadrant.
anglesL = (0.0, 45.0)
anglesR = (22.5, 67.5)
assert (all(0 <= x and x < 90 for x in anglesL + anglesR))
data = read_raw_data(raw_data_filename)
data = adjust_data_for_light_pulse_orientation(data)
dataL1R1, dataL1R2, dataL2R1, dataL2R2 = \
split_data_according_to_PBS_setting(data)
kappaL1R1 = \
compute_correlation_coefficient (anglesL[0], anglesR[0],
dataL1R1)
kappaL1R2 = \
compute_correlation_coefficient (anglesL[0], anglesR[1],
dataL1R2)
kappaL2R1 = \
compute_correlation_coefficient (anglesL[1], anglesR[0],
dataL2R1)
kappaL2R2 = \
compute_correlation_coefficient (anglesL[1], anglesR[1],
dataL2R2)
chsh_contrast = -kappaL1R1 + kappaL1R2 + kappaL2R1 + kappaL2R2
# The nominal value of the CHSH contrast for the chosen polarizer
# orientations is 2*sqrt(2).
chsh_contrast_nominal = 2 * sqrt(2.0)
print()
print(" light pulse events %9d" % len(data))
print()
print(" correlation coefs")
print(" %4.1f° %4.1f° %+9.6f" % (anglesL[0], anglesR[0],
kappaL1R1))
print(" %4.1f° %4.1f° %+9.6f" % (anglesL[0], anglesR[1],
kappaL1R2))
print(" %4.1f° %4.1f° %+9.6f" % (anglesL[1], anglesR[0],
kappaL2R1))
print(" %4.1f° %4.1f° %+9.6f" % (anglesL[1], anglesR[1],
kappaL2R2))
print()
print(" CHSH contrast %+9.6f" % chsh_contrast)
print(" 2*sqrt(2) = nominal %+9.6f" % chsh_contrast_nominal)
print(" difference %+9.6f" % (chsh_contrast -
chsh_contrast_nominal))
# A "CHSH violation" occurs if the CHSH contrast is >2.
# https://en.wikipedia.org/w/index.php?title=CHSH_inequality&oldid=1142431418
print()
print(" CHSH violation %+9.6f" % (chsh_contrast - 2))
print()
- Output:
An example using data from the Python implementation of Simulated optics experiment/Simulator.
light pulse events 100000
correlation coefs
0.0° 22.5° -0.707806
0.0° 67.5° +0.705415
45.0° 22.5° +0.712377
45.0° 67.5° +0.718882
CHSH contrast +2.844480
2*sqrt(2) = nominal +2.828427
difference +0.016053
CHSH violation +0.844480
RATFOR
There is a Ratfor preprocessor in C at https://sourceforge.net/p/chemoelectric/ratfor77 It outputs FORTRAN 77, which you can then compile like any FORTRAN 77 program.
Ratfor itself is a mix of Fortran syntax and C syntax. It was once, in a sense, the "Unix dialect of FORTRAN". (It was expanded into "Extended Fortran Language" but I do not remember what that added. I encountered EFL in System V Release 3.)
This program does the analysis in what might be the most efficient way possible. It uses very small constant space.
# -*- mode: indented-text; tab-width: 2; -*-
# Reference:
#
# A. F. Kracklauer, ‘EPR-B correlations: non-locality or geometry?’,
# J. Nonlinear Math. Phys. 11 (Supp.) 104–109 (2004).
# https://doi.org/10.2991/jnmp.2004.11.s1.13 (Open access, CC BY-NC)
# Standard FORTRAN 77 provides no way to allocate space at runtime.
# One would have to make a call to C, or do some such thing, or simply
# recompile if you need a bigger array. However, it is not
# necessary. This program computes the sums as it reads the data in.
# Once upon a time, because memory was precious and compilers made it
# hard to allocate memory, this would have been THE way to do such
# calculations. But the code ends up being less modular.
# In 2023 there is, for this particular program, no fathomable reason
# not to use double precision. So I shall use double precision.
define(angleL1, 0.0D0) # 'D' for double precision.
define(angleL2, 45.0D0)
define(angleR1, 22.5D0)
define(angleR2, 67.5D0)
subroutine count (N00, N01, N10, N11,
N00L1, N00L2, N00R1, N00R2,
N01L1, N01L2, N01R1, N01R2,
N10L1, N10L2, N10R1, N10R2,
N11L1, N11L2, N11R1, N11R2)
implicit none
integer N00, N01, N10, N11
integer N00L1, N00L2, N00R1, N00R2
integer N01L1, N01L2, N01R1, N01R2
integer N10L1, N10L2, N10R1, N10R2
integer N11L1, N11L2, N11R1, N11R2
integer fields(1:7)
integer events, i, tmp
continue
N00 = 0; N01 = 0; N10 = 0; N11 = 0
N00L1 = 0; N00L2 = 0; N00R1 = 0; N00R2 = 0
N01L1 = 0; N01L2 = 0; N01R1 = 0; N01R2 = 0
N10L1 = 0; N10L2 = 0; N10R1 = 0; N10R2 = 0
N11L1 = 0; N11L2 = 0; N11R1 = 0; N11R2 = 0
read (*,*) events
for (i = 0; i != events; i = i + 1)
{
read (*,*) fields
if (fields(1) == 1)
{
# Adjust for geometry.
tmp = fields(4); fields(4) = fields(5); fields(5) = tmp
tmp = fields(6); fields(6) = fields(7); fields(7) = tmp
}
if (fields(2) == 0 && fields(3) == 0)
{
N00 = N00 + 1
N00L1 = N00L1 + fields(4)
N00L2 = N00L2 + fields(5)
N00R1 = N00R1 + fields(6)
N00R2 = N00R2 + fields(7)
}
else if (fields(2) == 0 && fields(3) == 1)
{
N01 = N01 + 1
N01L1 = N01L1 + fields(4)
N01L2 = N01L2 + fields(5)
N01R1 = N01R1 + fields(6)
N01R2 = N01R2 + fields(7)
}
else if (fields(2) == 1 && fields(3) == 0)
{
N10 = N10 + 1
N10L1 = N10L1 + fields(4)
N10L2 = N10L2 + fields(5)
N10R1 = N10R1 + fields(6)
N10R2 = N10R2 + fields(7)
}
else
{
N11 = N11 + 1
N11L1 = N11L1 + fields(4)
N11L2 = N11L2 + fields(5)
N11R1 = N11R1 + fields(6)
N11R2 = N11R2 + fields(7)
}
}
end
function kappa (N, NL1, NL2, NR1, NR2)
implicit none
integer N, NL1, NL2, NR1, NR2
double precision kappa
double precision denom
double precision sinL, cosL, sinR, cosR
double precision sinLR, cosLR
continue
denom = N
# Equation (2.4) in the reference.
sinL = dsqrt (NL1 / denom)
cosL = dsqrt (NL2 / denom)
sinR = dsqrt (NR1 / denom)
cosR = dsqrt (NR2 / denom)
# Equation (2.3).
cosLR = (cosR * cosL) + (sinR * sinL)
sinLR = (sinR * cosL) - (cosR * sinL)
# Equation (2.5).
kappa = (cosLR * cosLR) - (sinLR * sinLR)
end
program OpSimA
implicit none
integer N00, N01, N10, N11
integer N00L1, N00L2, N00R1, N00R2
integer N01L1, N01L2, N01R1, N01R2
integer N10L1, N10L2, N10R1, N10R2
integer N11L1, N11L2, N11R1, N11R2
double precision kappa
double precision kappa00, kappa01, kappa10, kappa11
double precision chsh
continue
call count (N00, N01, N10, N11,
N00L1, N00L2, N00R1, N00R2,
N01L1, N01L2, N01R1, N01R2,
N10L1, N10L2, N10R1, N10R2,
N11L1, N11L2, N11R1, N11R2)
kappa00 = kappa (N00, N00L1, N00L2, N00R1, N00R2)
kappa01 = kappa (N01, N01L1, N01L2, N01R1, N01R2)
kappa10 = kappa (N10, N10L1, N10L2, N10R1, N10R2)
kappa11 = kappa (N11, N11L1, N11L2, N11R1, N11R2)
chsh = -kappa00 + kappa01 + kappa10 + kappa11
write (*,'()')
write (*,'(" light pulse events ", I9)') n00 + n01 + n10 + n11
write (*,'()')
write (*,'(" correlation coefs")')
write (*,'(" ", F4.1, "° ", F4.1, "° ", SPF9.6)') _
angleL1, angleR1, kappa00
write (*,'(" ", F4.1, "° ", F4.1, "° ", SPF9.6)') _
angleL1, angleR2, kappa01
write (*,'(" ", F4.1, "° ", F4.1, "° ", SPF9.6)') _
angleL2, angleR1, kappa10
write (*,'(" ", F4.1, "° ", F4.1, "° ", SPF9.6)') _
angleL2, angleR2, kappa11
write (*,'()')
write (*,'(" CHSH contrast ", SPF9.6)') chsh
write (*,'(" 2*sqrt(2) = nominal ", SPF9.6)') 2 * dsqrt (2.0D0)
write (*,'(" difference ", SPF9.6)') _
chsh - (2 * dsqrt (2.0D0))
write (*,'()')
write (*,'(" CHSH violation ", SPF9.6)') chsh - 2
write (*,'()')
end- Output:
Outputs may look different, due to implementation-dependent behavior of formatted output.
Here is output after compiling with gfortran:
light pulse events 100000
correlation coefs
0.0° 22.5° -0.707806
0.0° 67.5° +0.705415
45.0° 22.5° +0.712377
45.0° 67.5° +0.718882
CHSH contrast +2.844480
2*sqrt(2) = nominal +2.828427
difference +0.016053
CHSH violation +0.844480
But here is output after compiling with f2c:
light pulse events 100000
correlation coefs
.0° 22.5° -.707806
.0° 67.5° +.705415
45.0° 22.5° +.712377
45.0° 67.5° +.718882
CHSH contrast +2.844480
2*sqrt(2) = nominal +2.828427
difference +.016053
CHSH violation +.844480
Getting f2c formatted output to print what I would prefer to see is not worth my effort. Probably, if one is using f2c, it is better simply to remove the "SP" that add + signs ahead of numbers. Or, frankly, to hack libf2c more to one's liking. Maybe it was done this way simply to save printer ribbon ink.
Wren
Unlike the Python example, there is no option to read data from standard input.
import "io" for File
import "os" for Process
import "./seq" for Lst
import "./assert" for Assert
import "./fmt" for Fmt
/* The data for one light pulse. */
class PulseData {
construct new( logS, logL, logR, detectedL1, detectedL2, detectedR1, detectedR2) {
_logS = logS
_logL = logL
_logR = logR
_detectedL1 = detectedL1
_detectedL2 = detectedL2
_detectedR1 = detectedR1
_detectedR2 = detectedR2
}
logS { _logS }
logL { _logL }
logR { _logR }
detectedL1 { _detectedL1 }
detectedL2 { _detectedL2 }
detectedR1 { _detectedR1 }
detectedR2 { _detectedR2 }
toString {
return "PulseData(" + "%(_logS), %(_logL), %(_logR)" +
"%(_detectedL1), %(_detectedL2), %(_detectedR1), %(_detectedR2))"
}
// Swap detectedL1 and detectedL2.
swapLchannels() {
var t = _detectedL1
_detectedL1 = _detectedL2
_detectedL2 = t
}
// Swap detectedR1 and detectedR2.
swapRchannels() {
var t = _detectedR1
_detectedR1 = _detectedR2
_detectedR2 = t
}
// Swap channels on both right and left. This is done if the
// light source was (1,90°) instead of (1,0°). For in that case
// the orientations of the polarizing beam splitters, relative to
// the light source, is different by 90°.
swapLRchannels() {
swapLchannels()
swapRchannels()
}
}
// Split the data into two subsets, according to whether a set
// item satisfies a predicate. The return value is a tuple, with
// those items satisfying the predicate in the first tuple entry, the
// other items in the second entry.
var splitData = Fn.new { |predicate, data| Lst.partitions(data, predicate) }
// Some data items are for a (1,0°) light pulse. The others are
// for a (1,90°) light pulse. Thus the light pulses are oriented
// differently with respect to the polarizing beam splitters. We
// adjust for that distinction here.
var adjustDataForLightPulseOrientation = Fn.new { |data|
var pf = Fn.new { |item| item.logS == 0 }
var split = splitData.call(pf, data)
var data0 = split[0]
var data90 = split[1]
for (item in data90) item.swapLRchannels()
return data0 + data90
}
// Split the data into four subsets: one subset for each
// arrangement of the two polarizing beam splitters.
var splitDataAccordingToPBSSetting = Fn.new { |data|
var pfL = Fn.new { |item| item.logL == 0 }
var split = splitData.call(pfL, data)
var dataL1 = split[0]
var dataL2 = split[1]
var pfR = Fn.new { |item| item.logR == 0 }
var splitL1 = splitData.call(pfR, dataL1)
var dataL1R1 = splitL1[0]
var dataL1R2 = splitL1[1]
var splitL2 = splitData.call(pfR, dataL2)
var dataL2R1 = splitL2[0]
var dataL2R2 = splitL2[1]
return [dataL1R1, dataL1R2, dataL2R1, dataL2R2]
}
// Compute the correlation coefficient for the subset of the data
// that corresponding to a particular setting of the polarizing beam splitters.
var computeCorrelationCoefficient = Fn.new { |angleL, angleR, data|
// We make certain the orientations of beam splitters are
// represented by perpendicular angles in the first and fourth
// quadrant. This restriction causes no loss of generality, because
// the orientation of the beam splitter is actually a rotated "X".
Assert.ok([angleL, angleR].all { |x| 0 <= x && x < 90 })
var perpendicularL = angleL - 90 // In Quadrant 4.
var perpendicularR = angleR - 90 // In Quadrant 4.
// Note that the sine is non-negative in Quadrant 1, and the cosine
// is non-negative in Quadrant 4. Thus we can use the following
// estimates for cosine and sine. This is Equation (2.4) in the
// reference. (Note, one can also use Quadrants 1 and 2 and reverse
// the roles of cosine and sine. And so on like that.)
var N = data.count
var NL1 = 0
var NL2 = 0
var NR1 = 0
var NR2 = 0
for (item in data) {
NL1 = NL1 + item.detectedL1
NL2 = NL2 + item.detectedL2
NR1 = NR1 + item.detectedR1
NR2 = NR2 + item.detectedR2
}
var sinL = (NL1 / N).sqrt
var cosL = (NL2 / N).sqrt
var sinR = (NR1 / N).sqrt
var cosR = (NR2 / N).sqrt
// Now we can apply the reference's Equation (2.3).
var cosLR = (cosR * cosL) + (sinR * sinL)
var sinLR = (sinR * cosL) - (cosR * sinL)
// And then Equation (2.5).
return (cosLR * cosLR) - (sinLR * sinLR)
}
// Read the raw data into a list.
var readRawData = Fn.new { |filename|
var makeRecord = Fn.new { |line|
var x = line.split(" ")
x = x.map { |s| Num.fromString(s).truncate }.toList
return PulseData.new(x[0], x[1], x[2], x[3], x[4], x[5], x[6])
}
var data = []
var lines = File.read(filename).split("\n")
var numPulses = Num.fromString(lines[0])
for (i in 1..numPulses) {
data.add(makeRecord.call(lines[i]))
}
return data
}
var args = Process.arguments
if (args.count != 1) {
System.print("Please provide 1 command line argument: RAW_DATA_FILENAME")
return
}
var filename = args[0]
// Polarizing beam splitter orientations commonly used in actual
// experiments. These must match the values used in the simulation
// itself. They are by design all either zero degrees or in the
// first quadrant.
var anglesL = [0, 45]
var anglesR = [22.5, 67.5]
Assert.ok((anglesL + anglesR).all { |x| 0 <= x && x < 90 })
var data = readRawData.call(filename)
data = adjustDataForLightPulseOrientation.call(data)
var split = splitDataAccordingToPBSSetting.call(data)
var dataL1R1 = split[0]
var dataL1R2 = split[1]
var dataL2R1 = split[2]
var dataL2R2 = split[3]
var kappaL1R1 = computeCorrelationCoefficient.call(anglesL[0], anglesR[0], dataL1R1)
var kappaL1R2 = computeCorrelationCoefficient.call(anglesL[0], anglesR[1], dataL1R2)
var kappaL2R1 = computeCorrelationCoefficient.call(anglesL[1], anglesR[0], dataL2R1)
var kappaL2R2 = computeCorrelationCoefficient.call(anglesL[1], anglesR[1], dataL2R2)
var chshContrast = -kappaL1R1 + kappaL1R2 + kappaL2R1 + kappaL2R2
// The nominal value of the CHSH contrast for the chosen polarizer
// orientations is 2*sqrt(2).
var chshContrastNominal = 2.sqrt * 2
Fmt.print()
Fmt.print(" light pulse events $,10d", data.count)
Fmt.print()
Fmt.print(" correlation coefs")
Fmt.print(" $4.1f° $4.1f° $+9.6f", anglesL[0], anglesR[0], kappaL1R1)
Fmt.print(" $4.1f° $4.1f° $+9.6f", anglesL[0], anglesR[1], kappaL1R2)
Fmt.print(" $4.1f° $4.1f° $+9.6f", anglesL[1], anglesR[0], kappaL2R1)
Fmt.print(" $4.1f° $4.1f° $+9.6f", anglesL[1], anglesR[1], kappaL2R2)
Fmt.print()
Fmt.print(" CHSH contrast $+9.6f", chshContrast)
Fmt.print(" 2*sqrt(2) = nominal $+9.6f", chshContrastNominal)
Fmt.print(" difference $+9.6f", chshContrast - chshContrastNominal)
// A "CHSH violation" occurs if the CHSH contrast is > 2.
// https://en.wikipedia.org/w/index.php?title=CHSH_inequality&oldid=1142431418
Fmt.print()
Fmt.print(" CHSH violation $+9.6f", chshContrast - 2)
Fmt.print()
- Output:
An example using data from the Wren implementation of Simulated optics experiment/Simulator.
light pulse events 100,000
correlation coefs
0.0° 22.5° -0.702281
0.0° 67.5° +0.702892
45.0° 22.5° +0.696311
45.0° 67.5° +0.710296
CHSH contrast +2.811781
2*sqrt(2) = nominal +2.828427
difference -0.016646
CHSH violation +0.811781