jcs2opt.R

richard — Feb 21, 2014, 2:06 PM

## Port of Chantal Roth's simulation of Joy's model from Java to R

#### This variant: optimised for speed, now 10^6 runs
#### (optimisation is unfortunately at some cost to transparency)

## https://github.com/chenopodium/JCS
## https://github.com/chenopodium/JCS2

## See also Daniel Sabsay's Javascript version

## http://libertesphilosophica.info/eprsim/EPR_3-sphere_simulation5M.html
## http://libertesphilosophica.info/eprsim/eprsim.txt

## The formulas for Cab, Na, Nb are taken from Christian's 
## http://arxiv.org/abs/1301.1653
## Cab is given in (A27), Na and Nb are given by the expressions 
## just after (A28) and (A29).
## The "phases" phi_o^p etc. also come from this paper.
## The procedure generating an outcome +/-1 for the product of Alice and Bob's 
## measurement, or 0 for "no state", with the help of an auxiliary randomization, 
## is taken from Roth and Sabsay's code.

PHI_OP <- 0.000
PHI_OQ <- 0.000
PHI_OR <- -1.517
PHI_OS <- 0.663


set.seed(11091951)   ## set pseudo random generator's seed

angles <- seq(from = 0, to = 360, by = 7.5) * 2 * pi / 360

K <- length(angles)
corrs <- numeric(K)
Ns <- numeric(K)

M <- 10^6       ## number of runs

## generate M uniform random points on S^2
z <- runif(M, -1, 1)
t <- runif(M, 0, 2*pi)
r <- sqrt(1 - z^2)
x <- r * cos(t)
y <- r * sin(t)

es <- rbind(x, y, z)  ## 3 x M matrix
esUp <- rbind(z, x, y)
esDown <- rbind(y, z, x)

## generate M auxiliary uniforms for the biased coin tosses
unifs <- runif(M)

beta <- 0 * 2 * pi / 360  ## direction of measurement "b"

b <- c(cos(beta), sin(beta), 0)
bUp <- c(0, cos(beta), sin(beta))
bDown <- c(sin(beta), 0, cos(beta))



for (k in 1:(K-1)) {
    alpha <- angles[k]

    a <- c(cos(alpha), sin(alpha), 0)
    aUp <- c(0, cos(alpha), sin(alpha))
    aDown <- c(sin(alpha), 0, cos(alpha))

    S <- 0
    N <- 0

    eta_ae <- acos(colSums(es*a))
    eta_be <- acos(colSums(es*b))

    crossae <- esUp*aDown - esDown*aUp
    crossaeLengths <- sqrt(colSums(crossae^2))
    crossbe <- esUp*bDown - esDown*bUp
    crossbeLengths <- sqrt(colSums(crossbe^2))

    crossdot <- colSums(crossae*crossbe)/(crossaeLengths*crossbeLengths)

    N_a <- sqrt(cos(eta_ae + PHI_OP)^2 + sin(eta_ae + PHI_OQ)^2)
    N_b <- sqrt(cos(eta_be + PHI_OR)^2 + sin(eta_be + PHI_OS)^2)

    C_ab <- ( - cos(eta_ae + PHI_OP) * cos(eta_be + PHI_OR)
                + crossdot * sin(eta_ae + PHI_OQ) * sin(eta_be + PHI_OS)
            ) / (N_a * N_b)   

    good <- abs(C_ab) > unifs       ## Select the "states" 
    N <- sum(good) 
    corrs[k] <- sum(sign(C_ab[good]))/N
    Ns[k] <- N
}


corrs[K] <- corrs[1]
Ns[K] <- Ns[1]


plot(angles*180/pi, corrs, type="l", col="blue",
    main = "Two correlation functions",
    xlab = "Angle (degrees)",
    ylab = "Correlation")
points(angles*180/pi, corrs, col="blue", pch = ".", cex = 2)
lines(angles*180/pi, cos(angles), col = "black")
points(angles*180/pi, cos(angles), col = "black", pch = ".", cex = 2)

legend(x=0, y=-0.5, legend = c("Christian", "cosine"), text.col=c("blue", "black"), lty=1, col = c("blue", "black"))

plot(angles*180/pi, corrs, type="l", col="blue", xlim = c(0, 50), ylim = c(0.8, 1.0),
    main = "Two correlation functions",
    xlab = "Angle (degrees)",
    ylab = "Correlation")
points(angles*180/pi, corrs, type="b", col="blue", pch = ".", cex = 2)
lines(angles*180/pi, cos(angles), col = "black")
points(angles*180/pi, cos(angles), type = "b", col = "black", pch = ".", cex = 2)
legend(x = 0, y=0.85, legend = c("Christian", "cosine"), text.col=c("blue", "black"), lty=1, col = c("blue", "black"))

angles[1]*180/pi    # one particular angle

[1] 0

round(corrs[1], 4)  # simulation

[1] 0.9801

round(cos(angles[1]),  4)  # cosine

[1] 1

round(1/sqrt(Ns[1]), 4)  # standard error

[1] 0.0014