jcs.R

richard — Feb 20, 2014, 1:58 PM

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

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

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

angles <- seq(from = 0, to = 360, by = 15) * 2 * pi / 360 ## angles of measurement "a"

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

M <- 10^4       ## 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

## 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)

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


dot <- function(x, y) sum(x * y)       ## inner product


cross <- function (x,y) {            ## cross product
    c(x[2] * y[3] - x[3] * y[2],
        x[3] * y[1] - x[1] * y[3],
        x[1] * y[2] - x[2] * y[1]
    )
}

normalise <- function(x) {      ## normalise a non-zero vector to length 1
    l <- sqrt(sum(x^2))
    if (l > 0) x/l else x
}


for (k in 1:(K-1)) {
    alpha <- angles[k]
    a <- c(cos(alpha), sin(alpha), 0)
    S <- 0
    N <- 0

    for (run in 1:M) {
        e <- es[ , run]

        eta_ae <- acos(dot(a, e))
        eta_be <- acos(dot(b, e))

        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)

        crossae <- normalise(cross(a, e))    
        crossbe <- normalise(cross(b, e))

        C_ab <- ( - cos(eta_ae + PHI_OP) * cos(eta_be + PHI_OR)
                  + dot(crossae, crossbe) * sin(eta_ae + PHI_OQ) * sin(eta_be + PHI_OS)
                ) / (N_a * N_b)
        U <- unifs[run]
        if (U < abs(C_ab)) {
            S <- S + sign(C_ab)
            N <- N + 1
        }
    }
    sums[k] <- S
    Ns[k] <- N
}

sums[K] <- sums[1]
Ns[K] <- Ns[1]


plot(seq(from = 0, to = 360, by = 15), sums/Ns, ylim=c(-1, 1),
    main = "R port of Chantal Roth's JCS simulation",
    sub = "Number of runs = 10^4",
    xlab = "Angle",
    ylab = "Correlation"
)
lines(seq(from = 0, to = 360, by = 1), 
    cos(seq(from = 0, to = 360, by = 1) * 2 * pi / 360))

legend(0, -0.7, c("circles = simulation", "line = cosine"))