jcs2.R

richard — Feb 21, 2014, 12:47 PM

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

#### This variant: just one angle pair (alpha = beta = 0), but now 10^6 runs
#### TO DO: "vectorise" the algorithm to improve speed

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

K <- length(angles)
sums <- 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

## 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) {
    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/Ns ## correlation at alpha = beta = 0

[1] 0.9801