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.
 There was a typographical error in the formula for CHSH contrast, in which kappaL2R2 appeared where kappaL2R1 belonged. Both task description and Python code contained the error, which was discovered when Object Icon gave different results on the same raw data. Please review your examples to see if this error was copied. We apologize for the inconvenience.

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 ${\displaystyle 2{\sqrt {2}}\approx {2.828427}}$. 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.

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.

1. 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.
2. 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".
3. 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 ${\displaystyle \circ ~0~0~\circ ~\circ ~\circ ~\circ }$, one for lines that look like ${\displaystyle \circ ~0~1~\circ ~\circ ~\circ ~\circ }$, one for lines that look like ${\displaystyle \circ ~1~0~\circ ~\circ ~\circ ~\circ }$, and one for lines that look like ${\displaystyle \circ ~1~1~\circ ~\circ ~\circ ~\circ }$. Thus each group corresponds to a particular setting of the two polarizing beam splitters.
4. Separately compute a correlation coefficient for each group. This is done as follows.
1. Let ${\displaystyle N}$ equal the total number of lines in a group.
2. Let ${\displaystyle N_{L1}}$ 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.
3. Let ${\displaystyle N_{L2}}$ 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.
4. Add up similar numbers for the right side: ${\displaystyle N_{R1}}$ for the sixth entry in a line, and ${\displaystyle N_{R2}}$ for the seventh entry.
5. Compute the following numbers: ${\displaystyle sin_{L}={\sqrt {N_{L1}/N}}}$, ${\displaystyle cos_{L}={\sqrt {N_{L2}/N}}}$, ${\displaystyle sin_{R}={\sqrt {N_{R1}/N}}}$, and ${\displaystyle cos_{R}={\sqrt {N_{R2}/N}}}$. These numbers represent an estimate of the orientations of the two beam splitters, by inference from detection events in the light detectors.
6. Compute ${\displaystyle cos_{LR}=(cos_{R}\cdot cos_{L})+(sin_{R}\cdot sin_{L})}$ and ${\displaystyle sin_{LR}=(sin_{R}\cdot cos_{L})-(cos_{R}\cdot sin_{L})}$. These are angle difference formulas from trigonometry. We are "eliminating the coordinate system".
7. Compute ${\displaystyle \kappa =cos_{LR}^{2}-sin_{LR}^{2}}$. This is the correlation coefficient for the group, estimated by inference from detection events.
5. 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, and kappaL2R2.
6. Compute the estimated CHSH contrast: -kappaL1R1 + kappaL1R2 + kappaL2R1 + kappaL2R2.
7. 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 (${\displaystyle 2{\sqrt {2}}\approx {2.828427}}$) 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)

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)

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
{
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]. Translation of: Python 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 Translation of: Python 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:

if filename != "-":
let f = open(filename, fmRead)
f.close()
else:

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
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)
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
else
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")

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

--  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 Translation of: Python Library: Wren-seq Library: Wren-assert Library: Wren-fmt 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