## WARNING: Michel Fodje strongly disagrees my interpretation
## and use of the algorithms and mathematical model in his code.
## He has so far refused to comment on possible links with the
## "coincidence loophole" model Pascazio (1988), Larsson and Gill (2004),
## de Raedt, Hess, Michielsen (2011).
## The present code resulted from study of Michel Fodje's "epr-clocked"
## https://github.com/minkwe/epr-clocked (Python)
## Michel's algorithms are an instance of what has been known for quite
## a long time now as the "coincidence loophole".
## I have studied his code, extracted the mathematical algorithms
## and re-programmed part of them in R.
## This code attempts to reproduces the entire Python programme.
## See my "epr-clocked-core" for the code at the heart of the model.
## This version: June 2015. Original version: June 2014.
## Changes: *only* in this header text.
## Europhys. Lett., 67 (5), pp. 707–713 (2004)
## DOI: 10.1209/epl/i2004-10124-7
## Bell’s inequality and the coincidence-time loophole
## J.-A. Larsson and R. D. Gill
set.seed(-1) ## it took me a long time to find a seed
## which gave a violation of Larsson-Gill modified CHSH in one of these tests!
N <- 10^4 ## will make this larger, later
alpha <- c(0, 90) * pi / 180 # Alice's possible two settings
beta <- c(45, 135) * pi / 180 # Bob's possible two settings
a <- sample(c(1, 2), N, replace = TRUE) # Alice setting names (1, 2)
b <- sample(c(1, 2), N, replace = TRUE) # Bob setting names (1, 2)
coincWindow <- 0.0004
ts <- pi * 0.03 ## timescale
asym <- 0.98 ## asymmetry parameter
spin <- 0.5
n <- 2 * spin ## so actually, n = 1
phase <- pi * n ## so actually, phase = pi
el <- runif(N, 0, 2 * pi) ## minkwe's e (for left)
er <- el + phase ## minkwe's e + phase (for right)
p <- 0.5 * sin(runif(N, 0, pi / 6))^2
ml <- runif(N, asym, 1) ## small random jitter, left
mr <- runif(N, asym, 1) ## small random jitter, right
Cl <- (0.5/pi) * (-1)^n * cos(n * (el - alpha[a])) ## cos(e-a), left
Cr <- (0.5/pi) * (-1)^n * cos(n * (er - beta[b])) ## cos(e-a), right
tdl <- ts * pmax(ml * p - abs(Cl), 0) ## time delays, left
tdr <- ts * pmax(mr * p - abs(Cr), 0) ## time delays, right
A <- sign(Cl) ## measurement outcomes, left
B <- sign(Cr) ## measurement outcomes, right
AB <- A * B ## product of outcomes
mean(AB[a == 1 & b == 1 & abs(tdl-tdr) < coincWindow])
## [1] -0.6441785
mean(AB[a == 1 & b == 2 & abs(tdl-tdr) < coincWindow])
## [1] 0.6443954
mean(AB[a == 2 & b == 1 & abs(tdl-tdr) < coincWindow])
## [1] -0.6601307
mean(AB[a == 2 & b == 2 & abs(tdl-tdr) < coincWindow])
## [1] -0.6773109
S <- - mean(AB[a == 1 & b == 1 & abs(tdl-tdr) < coincWindow]) +
mean(AB[a == 1 & b == 2 & abs(tdl-tdr) < coincWindow]) -
mean(AB[a == 2 & b == 1 & abs(tdl-tdr) < coincWindow]) -
mean(AB[a == 2 & b == 2 & abs(tdl-tdr) < coincWindow])
S ## CHSH
## [1] 2.626015
## Now we compute the detection rate in order to check against the Larsson-Gill
## modified CHSH for coincidence selection statistics.
gamma11 <- sum((a == 1 & b == 1 & abs(tdl-tdr) < coincWindow)) /
sum((a == 1 & b == 1))
gamma12 <- sum((a == 1 & b == 2 & abs(tdl-tdr) < coincWindow)) /
sum((a == 1 & b == 2))
gamma21 <- sum((a == 2 & b == 1 & abs(tdl-tdr) < coincWindow)) /
sum((a == 2 & b == 1))
gamma22 <- sum((a == 2 & b == 2 & abs(tdl-tdr) < coincWindow)) /
sum((a == 2 & b == 2))
gamma <- min(c(gamma11, gamma12, gamma21, gamma22))
gamma
## [1] 0.7068459
2 + 6 *(1/gamma -1) ## Larsson-Gill bound (well: the bound is min of this and 4)
## [1] 4.488414
S ## Observed value of CHSH
## [1] 2.626015
## Now we superimpose the rest of Michel Fodje's model
## First, the emission times
EMIS_TIME_MAX <- 0.07 # maximum Time between emission events us
EMIS_TIME_MIN <- 0.014 # minimum Time between emission events us
EMIS_ALPHA <- 1.05
EMIS_BETA <- 6 - EMIS_ALPHA
emissionIntervals <- EMIS_TIME_MIN +
rbeta(N, EMIS_ALPHA, EMIS_BETA)*EMIS_TIME_MAX
emissionTimes <- cumsum(emissionIntervals)
measTimesLeft <- emissionTimes + tdl
measTimesRight <- emissionTimes + tdr
rc <- rnorm(N) ## this will determine lost particles
notLostLeft <- (rc < 4)
notLostRight <- (rc > -4)
TimesLeft <- measTimesLeft[notLostLeft]
OutcomesLeft <- A[notLostLeft]
SettingsLeft <- a[notLostLeft]
TimesRight <- measTimesRight[notLostRight]
OutcomesRight <- B[notLostRight]
SettingsRight <- b[notLostRight]
elementof <- function(el, set) {el %in% set}
elementofV <- Vectorize(elementof, "el")
Diffs <- abs(outer(TimesLeft, TimesRight, "-"))
diffs <- as.vector(Diffs)
LeftId <- row(Diffs)[Diffs < coincWindow]
RightId <- col(Diffs)[Diffs < coincWindow]
diffs <- Diffs[Diffs < coincWindow]
ordering <- order(diffs)
LeftPairs <- LeftId[ordering]
RightPairs <- RightId[ordering]
keep <- !(duplicated(LeftId) | duplicated(RightId))
LeftId <- LeftId[keep]
RightId <- RightId[keep]
Ao <- A[LeftPairs]
Bo <- B[RightPairs]
ao <- a[LeftPairs]
bo <- b[RightPairs]
ABo <- Ao * Bo ## product of outcomes
mean(ABo[ao == 1 & bo == 1])
## [1] -0.6441785
mean(ABo[ao == 1 & bo == 2])
## [1] 0.6443954
mean(ABo[ao == 2 & bo == 1])
## [1] -0.6601307
mean(ABo[ao == 2 & bo == 2])
## [1] -0.6773109
S <- - mean(ABo[ao == 1 & bo == 1]) +
mean(ABo[ao == 1 & bo == 2]) -
mean(ABo[ao == 2 & bo == 1]) -
mean(ABo[ao == 2 & bo == 2])
S
## [1] 2.626015
gammaL <- length(LeftPairs)/length(TimesLeft)
gammaR <- length(RightPairs)/length(TimesRight)
gamma <- pmin(gammaL, gammaR)
## Correct upper bound (Larsson-Gill, 2004)
2 + 6 * (1/gamma - 1)
## [1] 4.430518