Simulated optics experiment/Data analysis

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.
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

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 Object Icon and Python programs below can serve as a reference implementation.

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 , 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.
  4. Separately compute a correlation coefficient for each group. This is done as follows.
    1. Let equal the total number of lines in a group.
    2. 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.
    3. 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.
    4. Add up similar numbers for the right side: for the sixth entry in a line, and for the seventh entry.
    5. 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.
    6. Compute and . These are angle difference formulas from trigonometry. We are "eliminating the coordinate system".
    7. Compute . 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 () 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

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 + kappaL2R2 + 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

ObjectIcon

#!/bin/env 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"              %2.0r° %2.0r°   %+9.6r\n",
          angleL1, angleR1, kappaL1R1)
  printf (u"              %2.0r° %2.0r°   %+9.6r\n",
          angleL1, angleR2, kappaL1R2)
  printf (u"              %2.0r° %2.0r°   %+9.6r\n",
          angleL2, angleR1, kappaL2R1)
  printf (u"              %2.0r° %2.0r°   %+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° 23°   -0.699476
               0° 68°   +0.714172
              45° 23°   +0.692368
              45° 68°   +0.711478

        CHSH contrast   +2.817493
  2*sqrt(2) = nominal   +2.828427
           difference   -0.010935

       CHSH violation   +0.817493

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("              %2d° %2d°   %+9.6f" % (anglesL[0], anglesR[0],
                                                kappaL1R1))
    print("              %2d° %2d°   %+9.6f" % (anglesL[0], anglesR[1],
                                                kappaL1R2))
    print("              %2d° %2d°   %+9.6f" % (anglesL[1], anglesR[0],
                                                kappaL2R1))
    print("              %2d° %2d°   %+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° 22°   -0.707806
               0° 67°   +0.705415
              45° 22°   +0.712377
              45° 67°   +0.718882

        CHSH contrast   +2.844480
  2*sqrt(2) = nominal   +2.828427
           difference   +0.016053

       CHSH violation   +0.844480

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("              $2d° $2d°   $+9.6f", anglesL[0].floor, anglesR[0].floor, kappaL1R1)
Fmt.print("              $2d° $2d°   $+9.6f", anglesL[0].floor, anglesR[1].floor, kappaL1R2)
Fmt.print("              $2d° $2d°   $+9.6f", anglesL[1].floor, anglesR[0].floor, kappaL2R1)
Fmt.print("              $2d° $2d°   $+9.6f", anglesL[1].floor, anglesR[1].floor, 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° 22°   -0.702281
               0° 67°   +0.702892
              45° 22°   +0.696311
              45° 67°   +0.710296

        CHSH contrast   +2.811781
  2*sqrt(2) = nominal   +2.828427
           difference   -0.016646

       CHSH violation   +0.811781