A real EPR-Bohm scenario

[minkwe]

Conceptually, it is much better to lay out the mathematics of the model first (independent of any implementation) before trying to implement it in code, rather than try to figure out from a piece of code what the model is. That is why its been difficult for me to follow. So here is my summary of my understanding so far of the goal of the model in http://rpubs.com/jjc/233477:

Source:
- generates a random 3D vector e
- generates a random number f, such that f = -1 + 2/sqrt(1+3s/pi), where s is a uniform random number in the range [0, pi)

Station:
- z = +1 for Alice, -1 for Bob
- setting "a" (or "b") is a 3D vector constrained to the XY plane
- calculate the dot product "a.e"
- calculate q = {|a.e| if |a.e| > f else 0}
- detector output = z*sign("a.e")*sign(q)

Is that your understanding of the model? Forget for a moment about specific implementations and their verification.

[minkwe]

So I am still trying to decipher if your Python code is actually verifying the R code. You never answered my question about returning "z". I still don't understand that.

First let us agree on what the essential mathematics of the model is, then we can ask the question whether the R/Python/Mathematica code correctly implements it or not. I don't understand how to "verify the R code" any other way. The model as expressed in my previous post is completely agnostic to the way the settings "a" and "b" are chosen. It could be random or sequential and it won't make a difference. If it appears to, then there is an error in the code.

[Joy Christian]

So I am still trying to decipher if your Python code is actually verifying the R code. You never answered my question about returning "z". I still don't understand that.

First let us agree on what the essential mathematics of the model is, then we can ask the question whether the R/Python/Mathematica code correctly implements it or not. I don't understand how to "verify the R code" any other way. The model as expressed in my previous post is completely agnostic to the way the settings "a" and "b" are chosen. It could be random or sequential and it won't make a difference. If it appears to, then there is an error in the code.

In this R code the function g = the function p x the function q. So it seems to me to be completely equivalent to this R code.

So how can there be an error? At least mathematically there can't be an error if we write g = p x q and use p x q instead of g.

PS: Note also that sign(X x Y) = sign(X) x sign(Y).

***

[FrediFizzx]

Conceptually, it is much better to lay out the mathematics of the model first (independent of any implementation) before trying to implement it in code, rather than try to figure out from a piece of code what the model is. That is why its been difficult for me to follow. So here is my summary of my understanding so far of the goal of the model in http://rpubs.com/jjc/233477:

Source:
- generates a random 3D vector e
- generates a random number f, such that f = -1 + 2/sqrt(1+3s/pi), where s is a uniform random number in the range [0, pi)

Station:
- z = +1 for Alice, -1 for Bob
- setting "a" (or "b") is a 3D vector constrained to the XY plane
- calculate the dot product "a.e"
- calculate q = {|a.e| if |a.e| > f else 0}
- detector output = z*sign("a.e")*sign(q)

Is that your understanding of the model? Forget for a moment about specific implementations and their verification.

That is right if you take out the "z" part and change the first |a.e| to cos(a) in the q-function.

We have separated out the polarizer to sign(a.e) so all we need for the q-function is cos(a). Now, the model that works in Mathematica is similar but with the other polarizer action of |a.e|<pi returns +1 else returns -1. And that also works in R. So we have,

Station:
- setting "a" (or "b") is a 3D vector constrained to the XY plane
- calculate the dot product "a.e" = c
- calculate polarizer = {+1 if |c|<pi else -1}
- calculate q = {cos(a*polarizer) if |c| > f else 0}
- detector output = sign(q)

[minkwe]

Station:
- setting "a" (or "b") is a 3D vector constrained to the XY plane
- calculate the dot product "a.e" = c
- calculate polarizer = {+1 if |c|<pi else -1}
- calculate q = {cos(a*polarizer) if |c| > f else 0}
- detector output = sign(q)

This makes no sense to me. What does "polarizer" mean physically? What do you think "a" and "b" represent? "c" is a cosine which means it ranges between -1 and +1. It will always be less than pi, so your "polarizer" will always be +1. In any case, the above is NOT what the code http://rpubs.com/jjc/233477 is doing. The cos(Angles) in the code is problematic. It can't be right.

[minkwe]

So I am still trying to decipher if your Python code is actually verifying the R code. You never answered my question about returning "z". I still don't understand that.

First let us agree on what the essential mathematics of the model is, then we can ask the question whether the R/Python/Mathematica code correctly implements it or not. I don't understand how to "verify the R code" any other way. The model as expressed in my previous post is completely agnostic to the way the settings "a" and "b" are chosen. It could be random or sequential and it won't make a difference. If it appears to, then there is an error in the code.

In this R code the function g = the function p x the function q. So it seems to me to be completely equivalent to this R code.

So how can there be an error? At least mathematically there can't be an error if we write g = p x q and use p x q instead of g.

***

Joy, the devil is in the details.

Code: Select all

p = function(u,v){colSums(u*v)} # polarizer function of Bell
q = function(u,v,s){ifelse(abs(p(u,v)) > f, cos(Angles), 0)}
g = function(u,v,s){p(u,v)*q(u,v,s)}   # cf. http://rpubs.com/jjc/84238

vs

Code: Select all

g = function(u,v,s){ifelse(abs(colSums(u*v)) > f, colSums(u*v), 0)}

Where is the "Angles" in "g"?


My understanding is that the original model was:

Station:
- z = +1 for Alice, -1 for Bob
- setting "a" (or "b") is a 3D vector constrained to the XY plane
- calculate the dot product "a.e"
- detector output = z*sign("a.e") if |a.e| > f else 0

Isn't the goal is to add "polarizer action" to this without changing the essential core of the model? I dont think the model such as

Station:
- setting "a" (or "b") is a 3D vector constrained to the XY plane
- calculate the dot product "a.e" = c
- calculate polarizer = {+1 if |c|<pi else -1}
- calculate q = {cos(a*polarizer) if |c| > f else 0}
- detector output = sign(q)

is doing that. It is completely different.

Fred,
z is just a generalization of the fact that we have to multiply by -1 on Bob's side. Rather than write two separate functions with one multiplied by -1, we simply define z = (+1 for Alice, -1 for Bob) and write down a single function.

[FrediFizzx]

Station:
- setting "a" (or "b") is a 3D vector constrained to the XY plane
- calculate the dot product "a.e" = c
- calculate polarizer = {+1 if |c|<pi else -1}
- calculate q = {cos(a*polarizer) if |c| > f else 0}
- detector output = sign(q)

This makes no sense to me. What does "polarizer" mean physically? What do you think "a" and "b" represent? "c" is a cosine which means it ranges between -1 and +1. It will always be less than pi, so your "polarizer" will always be +1. In any case, the above is NOT what the code http://rpubs.com/jjc/233477 is doing. The cos(Angles) in the code is problematic. It can't be right.

Oops, sorry. Right, c = |a - e| not cos(a - e). Let me try again.

Station:
- setting "a" (or "b") is a 3D vector constrained to the XY plane
- calculate |a - e| = c
- calculate polarizer = {+1 if c<pi else -1}
- calculate q = {cos(a*polarizer) if c > f else 0}
- detector output = sign(q)[/quote]

IOW, we are just using the angle part of the vectors.
.

[minkwe]

Oops, sorry. Right, c = |a - e| not cos(a - e). Let me try again.

Station:
- setting "a" (or "b") is a 3D vector constrained to the XY plane
- calculate |a - e| = c
- calculate polarizer = {+1 if c<pi else -1}
- calculate q = {cos(a*polarizer) if c > f else 0}
- detector output = sign(q)

IOW, we are just using the angle part of the vectors.
.

This is a completely different model from the original one and not even the same as what Joy is doing in the new R-code. Not to mention that you are mixing 2D and 3D representations. "a" and "b" and "e" are vectors not angles.

[minkwe]

Let me stick to 2D to see if I understand what this "new" model is:

Station:
- setting "a" (or "b") is a random angle between (0, 2pi), e is also an angle between (0, 2pi)
- calculate |a - e| = c
- calculate polarizer = {+1 if c<pi else -1}
- calculate q = {cos(a*polarizer) if c > f else 0}
- detector output = sign(q)

Is that the model?

[FrediFizzx]

Oops, sorry. Right, c = |a - e| not cos(a - e). Let me try again.

Station:
- setting "a" (or "b") is a 3D vector constrained to the XY plane
- calculate |a - e| = c
- calculate polarizer = {+1 if c<pi else -1}
- calculate q = {cos(a*polarizer) if c > f else 0}
- detector output = sign(q)

IOW, we are just using the angle part of the vectors.
.

This is a completely different model from the original one and not even the same as what Joy is doing in the new R-code. Not to mention that you are mixing 2D and 3D representations. "a" and "b" and "e" are vectors not angles.

Ok right, the Mathematica version uses angles and I forgot to mention that "e" is 2D in it. Let me find the R version that implements "calculate polarizer = {+1 if c<pi else -1}".

[Joy Christian]

Joy, the devil is in the details.

Code: Select all

p = function(u,v){colSums(u*v)} # polarizer function of Bell
q = function(u,v,s){ifelse(abs(p(u,v)) > f, cos(Angles), 0)}
g = function(u,v,s){p(u,v)*q(u,v,s)}   # cf. http://rpubs.com/jjc/84238

vs

Code: Select all

g = function(u,v,s){ifelse(abs(colSums(u*v)) > f, colSums(u*v), 0)}

Where is the "Angles" in "g"?

If I use

Code: Select all

g = function(u,v,s){p(u,v)*q(u,v,s)}

instead of

Code: Select all

g = function(u,v,s){ifelse(abs(colSums(u*v)) > f, colSums(u*v), 0)}

in the original R code, then it runs just fine and produces the strong correlations.

***

[FrediFizzx]

OK, here is the R code that works with the same polarizer function as Mathematica.

Code: Select all

Angles = seq(from = 0, to = 360, by = 7.2) * 2 * pi/360

K = length(Angles) # The total number of angles between 0 and 2pi

corrs = matrix(nrow = K, ncol = K, data = 0) # Container for correlations

M = 10^4 # Size of the pre-ensemble. Next one can try 10^5, or even 10^6

r = runif(M, 0, 2*pi) # M uniformly distributed numbers between 0 and 2pi

z = runif(M, -1, +1) # M uniformly distributed numbers between -1 and +1

h = sqrt(1 - z^2)

x = h * cos(r)

y = h * sin(r)

e1 = rbind(x, y, z)
e2 = -e1

s = runif(M, 0, pi) # Initial states of the spins are the pairs (e, s) within S^3

f = -1 + (2/sqrt(1 + ((3 * s)/pi))) # For details see the paper arXiv:1405.2355

g1 = function(u1,v1){ifelse(abs(colSums(u1*v1)) > f, cos(Angles), 0)}
g2 = function(u2,v2){ifelse(abs(colSums(u2*v2)) > f, cos(Angles), 0)}

d = function(o,t){ifelse(2*acos(colSums(o*t)) < pi, 1, -1)}

for (i in 1:K) {
  # i = 5
  alpha = Angles[i]
  a = c(cos(alpha), sin(alpha), 0)  # Measurement direction 'a'
 
  for (j in 1:K) {
    # j = 23   
    beta = Angles[j]
    b = c(cos(beta), sin(beta), 0)  # Measurement direction 'b'
    A = sign(d(a,e1)*(g1(a,e1)))  # Alice's measurement results A(a, e, s) = +/-1
    B = sign(d(b,e2)*(g2(b,e2)))  # Bob's measurement results B(b, e, s) = -/+1
   
    N = length((A*B)[A & B]) # Number of all possible events observed in S^3
   
    corrs[i,j] = sum(A*B)/N  # Product moment correlation coefficient E(a, b)
   
  }
}

par(mar = c(0, 0, 2, 0))
persp(x = Angles, y = Angles, main = "The strong correlations predicted by the 3-sphere model", z = corrs, zlim = c(-1, 1), col = "pink", theta = 135, phi = 30, scale = FALSE, xlab = "alpha", ylab = "beta")

This following line is the polarizer function.

d = function(o,t){ifelse(2*acos(colSums(o*t)) < pi, 1, -1)}
.

[minkwe]

Let me stick to 2D to see if I understand what this "new" model is:

Station:
- setting "a" (or "b") is a random angle between (0, 2pi), e is also an angle between (0, 2pi)
- calculate |a - e| = c
- calculate polarizer = {+1 if c<pi else -1}
- calculate q = {cos(a*polarizer) if c > f else 0}
- detector output = sign(q)

Is that the model?

Fred could you confirm if the above is correct in 2D? In effect the polarizer calculation is

sign(pi-c)

If the above is correct, then I can do the 2D test of the exact model above in python. I prefer to go from a clear agreement of what the model is, rather than trying to read the R-code.

[minkwe]

If I use

Code: Select all

g = function(u,v,s){p(u,v)*q(u,v,s)}

instead of

Code: Select all

g = function(u,v,s){ifelse(abs(colSums(u*v)) > f, colSums(u*v), 0)}

in the original R code, then it runs just fine and produces the strong correlations.

***

The problem is with the q function. What is the meaning of "Angles" in that expression

Code: Select all

q = function(u,v,s){ifelse(abs(p(u,v)) > f, cos(Angles), 0)}

u and v are unit vectors that change from one iteration to the next. s is a scalar, f is also a scalar. p is the dot product which is also a scalar, Angles is a vector of dimention 51 as defined earlier on the first line. What does cos(Angles) mean then??

[minkwe]

Let me stick to 2D to see if I understand what this "new" model is:

Station:
- setting "a" (or "b") is a random angle between (0, 2pi), e is also an angle between (0, 2pi)
- calculate |a - e| = c
- calculate polarizer = {+1 if c<pi else -1}
- calculate q = {cos(a*polarizer) if c > f else 0}
- detector output = sign(q)

Is that the model?

Fred could you confirm if the above is correct in 2D? In effect the polarizer calculation is

sign(pi-c)

If the above is correct, then I can do the 2D test of the exact model above in python. I prefer to go from a clear agreement of what the model is, rather than trying to read the R-code.

If you agree that the above is the model, then it does not work. Here is what the python simulation produces for the above 2D.

https://drive.google.com/open?id=0B6sZy ... EdWcTFoNUE

Code: Select all

def particle_emitter():
    e = numpy.random.uniform(0, 2*numpy.pi)
    s = numpy.random.uniform(1, 4)
    f = -1 + 2 / numpy.sqrt(s)
    Pr = (e, f, 1)  # Right particle
    Pl = (e, f, -1)  # Left particle
    return Pl, Pr


def particle_detector(info):
    """Calculate and return the station outcome for the given `particle` and setting"""
    particle, setting = info
    a = setting
    e, f, z = particle
    c = abs(a-e)
    pol = numpy.sign(numpy.pi - c)
    q   = numpy.cos(a*pol) if c > f else 0
    out = z*numpy.sign(q)

    return [0, setting,out]


[Joy Christian]

If I use

Code: Select all

g = function(u,v,s){p(u,v)*q(u,v,s)}

instead of

Code: Select all

g = function(u,v,s){ifelse(abs(colSums(u*v)) > f, colSums(u*v), 0)}

in the original R code, then it runs just fine and produces the strong correlations.

***

The problem is with the q function. What is the meaning of "Angles" in that expression

Code: Select all

q = function(u,v,s){ifelse(abs(p(u,v)) > f, cos(Angles), 0)}

u and v are unit vectors that change from one iteration to the next. s is a scalar, f is also a scalar. p is the dot product which is also a scalar, Angles is a vector of dimention 51 as defined earlier on the first line. What does cos(Angles) mean then??

I checked cos(Angles) in the code and got the following answer:

Code: Select all

(cos(Angles))
##  [1]  1.00000000  0.99211470  0.96858316  0.92977649  0.87630668
##  [6]  0.80901699  0.72896863  0.63742399  0.53582679  0.42577929
## [11]  0.30901699  0.18738131  0.06279052 -0.06279052 -0.18738131
## [16] -0.30901699 -0.42577929 -0.53582679 -0.63742399 -0.72896863
## [21] -0.80901699 -0.87630668 -0.92977649 -0.96858316 -0.99211470
## [26] -1.00000000 -0.99211470 -0.96858316 -0.92977649 -0.87630668
## [31] -0.80901699 -0.72896863 -0.63742399 -0.53582679 -0.42577929
## [36] -0.30901699 -0.18738131 -0.06279052  0.06279052  0.18738131
## [41]  0.30901699  0.42577929  0.53582679  0.63742399  0.72896863
## [46]  0.80901699  0.87630668  0.92977649  0.96858316  0.99211470
## [51]  1.00000000


[minkwe]

As I expected, so for each iteration, the q function is returning the exact same thing when p > f. It is returning 51 numbers. So g= p*q will also be 51 numbers in those cases and also sign(q) will be 51 numbers some -1 other +1 and so will sign(p)*sign(q) etc. Do you see now that this cannot be correct?

[FrediFizzx]

As I expected, so for each iteration, the q function is returning the exact same thing when p > f. It is returning 51 numbers. So g= p*q will also be 51 numbers in those cases and also sign(q) will be 51 numbers some -1 other +1 and so will sign(p)*sign(q) etc. Do you see now that this cannot be correct?

No it is not returning 51 numbers. It is 51 x 51 because the iterations are nested.

[Joy Christian]

As I expected, so for each iteration, the q function is returning the exact same thing when p > f. It is returning 51 numbers. So g= p*q will also be 51 numbers in those cases and also sign(q) will be 51 numbers some -1 other +1 and so will sign(p)*sign(q) etc. Do you see now that this cannot be correct?

Not really. But what is more confusing is: why does the code work at all if q is wrong? If I remove q by setting f = 0, then I get the linear correlations. So somehow the q-factor produces the right correlation. Don't you find that curious? If q is wrong (and I am sure you are right about that), then why does p x q work in place of original g, with linear background restored as expected when q is set to 1 (by setting f =0)?

***

[minkwe]

As I expected, so for each iteration, the q function is returning the exact same thing when p > f. It is returning 51 numbers. So g= p*q will also be 51 numbers in those cases and also sign(q) will be 51 numbers some -1 other +1 and so will sign(p)*sign(q) etc. Do you see now that this cannot be correct?

No it is not returning 51 numbers. It is 51 x 51 because the iterations are nested.

Each time the function is called it returns 51 numbers. The R code calls it 51x51 times. You can add a print statement inside the inner loop to print

length(q(a,e,s)) or length(sign(q(a,e,s)))

and you will see that it gives 51 each time for a total of 51x51 times