<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Introduction
This notebook is the computational appendix of arXiv
Step1:  Then we define the scenario we are considering. In full generality we define it with the three parameters $(N,m,d)$, corresponding to the case of $N$ parties, $m$ measurement choices with $d$ outcomes each. For instance
Step2:  We generated the symbolic variables defining the measurements. They are treated as commuting variables, since to test for locality, we use the NPA hierarchy with commuting measurements. Given that we will always work in the correlator space, we define the substitution rule $\mathcal({M}_{x_i}^{(i)})^2 = \mathbb{1}$ for any $i = 1,...,N$ and $x_i = 0,..,m-1$.
Step3:  Feasibility problem
First, we introduce the code to solve the problem of determinining whether a set of correlation is nonlocal. We cast it as an SDP feasibility problem, so that whenever one gets an infeasible solution it implies that the given correlations cannot be described by a local model.
$W$ state example
We start with generating the state and operators corresponding the measurements
Step4:  Then we generate the moments, i.e. the values of the correlators
Step5:  We construct the hierarchy at level two with corresponding substitutions and moments assignments
Step6:  Finally, we solve the SDP feasibility problem
Step7:  $GHZ$ state example
We start with generating the state and operators corresponding the measurements
Step8:  For the $GHZ$ state, we need to add the values of the two full-body correlators $\langle M_0^{(1)}M_0^{(2)}...M_0^{(N)} \rangle$ and $\langle M_1^{(1)}M_0^{(2)}...M_0^{(N)} \rangle$.
Step9:  We construct the hierarchy at the hybrid level $\lbrace \mathcal{O}_2,M_0^{(1)}M_0^{(2)}...M_0^{(N)},M_1^{(1)}M_0^{(2)}...M_0^{(N)} \rbrace$ with corresponding substitutions and moments assignments
Step10:  We solve the SDP feasibility problem
Step11:  Noise robustness
We show how to estimate the noise robustness of some observed nonlocal correlations numerically. Given a set of correlations generated by a state $\rho$ and some measurements $M_{x_i}^{(i)}$, we estimate the minimal amount of noise $\lambda$ such that the correlations arise from $(1 - \lambda )\rho + \lambda \frac{\mathbb{1}}{2^N}$.
If we perform measurements on a Bloch sphere, the noise affects the values of the correlators the same way. That is, if we denote $v_{j_1 j_2 j_3 j_4}^{i_1 i_2 i_3 i_4}$ the value of $\langle M_{j_1}^{(i_1)}M_{j_2}^{(i_2)}M_{j_3}^{(i_3)}M_{j_4}^{(i_4)} \rangle$, the noisy state will result in the corresponding value of
$(1 - \lambda ) v_{j_1 j_2 j_3 j_4}^{i_1 i_2 i_3 i_4} $
We cast this problem as an SDP by introducing $\lambda$ as the objective function of the minimization problem and by replacing the moments substitutions with the modified values depending on $\lambda$. First, we define $\lambda$ as an additional symbolic variable
Step12:  $W$ state example
We generate the moment constraints as a function of the $\lambda$ parameter
Step13:  Then the relaxation
Step14:  and finally we solve the SDP
Step15:  $GHZ$ state example
We repeat the same procedure as for the $W$ state case, with the addition of the two full-body correlators
Step16:  The relaxation is generated as
Step17:  Solving it, we get
Step18:  Dual inequality
From the dual of the noise robustness SDP, one can extract a Bell inequality that detects the observed nonlocal correlations. The derived inequality will be in the form $\sum \alpha_{j_1 j_2 j_3 j_4}^{i_1 i_2 i_3 i_4} \langle M_{j_1}^{(i_1)}M_{j_2}^{(i_2)}M_{j_3}^{(i_3)}M_{j_4}^{(i_4)} \rangle + 1 \geq 0$.
We will present the code we used for the case of the $GHZ$ and $m = 2$ measurement choices. First, we solve the noise robustness SDP
Step19:  We only assign the values of up to the two-body correlators. This will ensure that we get the same expression for the Bell inequality as the one presented in the manuscript. 
Now, we proceed to generate a second relaxation, but we do not substitute the monomials corresponding to the observed correlations.
Step20:  Then we assign to the second SDP the same dual variables as for the one solved before.
Step21:  We extract the value of the dual variable for each monomial appearing in the moment matrix. For the function "extract_dual_value" to work properly we needed to generate a moment matrix in which the values of the correlators where not substituted yet.
Step22:  Lastly, we normalize the inequality so to be in the form presented before (i.e. with the classical bound being equal to $1$ )
<ASSISTANT_TASK:>
Python Code:
from __future__ import print_function, division
from math import sqrt
from qutip import sigmax, sigmaz
from ncpol2sdpa import flatten, SdpRelaxation, generate_variables
from time import time
from sympy import S
from local_tools import generate_commuting_measurements, get_W_reduced, \
    get_GHZ_reduced, get_moment_constraints, get_fullmeasurement
Explanation: Introduction
This notebook is the computational appendix of arXiv:1612.08551. We demonstrate how to use a modified version of the NPA hierarchy to test whether a set of observed correlations is local. We also give the details how to reproduce the numerical results shown in the paper. Furthermore, we show how to determine the robustness to noise of a set of given nonlocal correlations and to extract a Bell inequality from the dual of the SDP.
To improve readability of this notebook, we placed the supporting functions to a separate file; please download this in the same folder as the notebook if you would like to evaluate it. The following dependencies must also be available: at least one SDP solver (SDPA as an executable in the path or Mosek with its Python interface installed; cvxopt as a solver is not recommended) together with the Ncpol2sdpa and qutip packages.
First, we import everything we will need:
End of explanation
N, m, d = 7, 2, 2
configuration = [d for _ in range(m)]
measurements = [generate_commuting_measurements(configuration, chr(65+i))
                for i in range(N)]
Explanation: Then we define the scenario we are considering. In full generality we define it with the three parameters $(N,m,d)$, corresponding to the case of $N$ parties, $m$ measurement choices with $d$ outcomes each. For instance
End of explanation
substitutions = {M**2: S.One for M in flatten(measurements)}
Explanation: We generated the symbolic variables defining the measurements. They are treated as commuting variables, since to test for locality, we use the NPA hierarchy with commuting measurements. Given that we will always work in the correlator space, we define the substitution rule $\mathcal({M}_{x_i}^{(i)})^2 = \mathbb{1}$ for any $i = 1,...,N$ and $x_i = 0,..,m-1$.
End of explanation
W_state = get_W_reduced(N)
W_operators = [[sigmax(), sigmaz()] for _ in range(N)]
Explanation: Feasibility problem
First, we introduce the code to solve the problem of determinining whether a set of correlation is nonlocal. We cast it as an SDP feasibility problem, so that whenever one gets an infeasible solution it implies that the given correlations cannot be described by a local model.
$W$ state example
We start with generating the state and operators corresponding the measurements:
End of explanation
time0 = time()
moments = get_moment_constraints(N, W_state, measurements, W_operators)
print("Constraints were generated in " + str(time()-time0) + " seconds.")
Explanation: Then we generate the moments, i.e. the values of the correlators:
End of explanation
time0 = time()
sdp = SdpRelaxation(flatten(measurements))
sdp.get_relaxation(2, substitutions=substitutions, momentsubstitutions=moments)
print("SDP relaxation was generated in " + str(time()-time0) + " seconds.")
Explanation: We construct the hierarchy at level two with corresponding substitutions and moments assignments:
End of explanation
sdp.solve(solver="mosek")
print("SDP was solved in " + str(sdp.solution_time) + " seconds.")
print("SDP status is " + sdp.status)
Explanation: Finally, we solve the SDP feasibility problem:
End of explanation
GHZ_state = get_GHZ_reduced(4)
GHZ_operators = [[sigmax(), 1/sqrt(2)*(sigmaz()+sigmax())] for _ in range(N)]
Explanation: $GHZ$ state example
We start with generating the state and operators corresponding the measurements:
End of explanation
time0 = time()
moments = get_moment_constraints(N, GHZ_state, measurements, GHZ_operators)
extra = get_fullmeasurement([0 for _ in range(N)], measurements)
extramonomials = [extra]
moments[extra] = 1
extra = get_fullmeasurement(flatten([1,[0 for _ in range(N-1)]]), measurements)
extramonomials.append(extra)
moments[extra] = (1/sqrt(2))
print("Constraints were generated in " + str(time()-time0) + " seconds.")
Explanation: For the $GHZ$ state, we need to add the values of the two full-body correlators $\langle M_0^{(1)}M_0^{(2)}...M_0^{(N)} \rangle$ and $\langle M_1^{(1)}M_0^{(2)}...M_0^{(N)} \rangle$.
End of explanation
time0 = time()
sdp = SdpRelaxation(flatten(measurements), verbose=1, parallel=True)
sdp.get_relaxation(2, substitutions=substitutions, 
                   momentsubstitutions=moments, extramonomials=extramonomials)
print("SDP relaxation was generated in " + str(time()-time0) + " seconds.")
Explanation: We construct the hierarchy at the hybrid level $\lbrace \mathcal{O}_2,M_0^{(1)}M_0^{(2)}...M_0^{(N)},M_1^{(1)}M_0^{(2)}...M_0^{(N)} \rbrace$ with corresponding substitutions and moments assignments:
End of explanation
sdp.solve(solver="mosek")
print("SDP was solved in " + str(sdp.solution_time) + " seconds.")
print("SDP status is " + sdp.status)
Explanation: We solve the SDP feasibility problem
End of explanation
lambda_ = generate_variables("\lambda")[0]
Explanation: Noise robustness
We show how to estimate the noise robustness of some observed nonlocal correlations numerically. Given a set of correlations generated by a state $\rho$ and some measurements $M_{x_i}^{(i)}$, we estimate the minimal amount of noise $\lambda$ such that the correlations arise from $(1 - \lambda )\rho + \lambda \frac{\mathbb{1}}{2^N}$.
If we perform measurements on a Bloch sphere, the noise affects the values of the correlators the same way. That is, if we denote $v_{j_1 j_2 j_3 j_4}^{i_1 i_2 i_3 i_4}$ the value of $\langle M_{j_1}^{(i_1)}M_{j_2}^{(i_2)}M_{j_3}^{(i_3)}M_{j_4}^{(i_4)} \rangle$, the noisy state will result in the corresponding value of
$(1 - \lambda ) v_{j_1 j_2 j_3 j_4}^{i_1 i_2 i_3 i_4} $
We cast this problem as an SDP by introducing $\lambda$ as the objective function of the minimization problem and by replacing the moments substitutions with the modified values depending on $\lambda$. First, we define $\lambda$ as an additional symbolic variable:
End of explanation
time0 = time()
moments = get_moment_constraints(N, W_state, measurements, W_operators, lambda_)
print("Constraints were generated in " + str(time()-time0) + " seconds.")
Explanation: $W$ state example
We generate the moment constraints as a function of the $\lambda$ parameter:
End of explanation
sdp = SdpRelaxation(flatten(measurements), parameters=[lambda_],
                    verbose=1, parallel=True)
sdp.get_relaxation(2, objective=lambda_, substitutions=substitutions, 
                   momentsubstitutions=moments)
print("SDP relaxation was generated in " + str(time()-time0) + " seconds.")
Explanation: Then the relaxation:
End of explanation
sdp.solve(solver="mosek")
print("SDP was solved in " + str(sdp.solution_time) + " seconds.")
print("lambda_min is " + str(sdp.primal))
Explanation: and finally we solve the SDP:
End of explanation
time0 = time()
moments = get_moment_constraints(N, GHZ_state, measurements, GHZ_operators, lambda_)
extra = get_fullmeasurement([0 for _ in range(N)], measurements)
extramonomials = [extra]
moments[extra] = (1 - lambda_)
extra = get_fullmeasurement(flatten([1, [0 for _ in range(N-1)]]), measurements)
extramonomials.append(extra)
moments[extra] = (1/sqrt(2))*(1 - lambda_)
print("Constraints were generated in " + str(time()-time0) + " seconds.")
Explanation: $GHZ$ state example
We repeat the same procedure as for the $W$ state case, with the addition of the two full-body correlators
End of explanation
sdp = SdpRelaxation(flatten(measurements), parameters=[lambda_],verbose=1, parallel=True)
sdp.get_relaxation(2, objective=lambda_, substitutions=substitutions, 
                   momentsubstitutions=moments, extramonomials=extramonomials)
print("SDP relaxation was generated in " + str(time()-time0) + " seconds.")
Explanation: The relaxation is generated as:
End of explanation
sdp.solve(solver="mosek")
print("SDP was solved in " + str(sdp.solution_time) + " seconds.")
print("lambda_min is " + str(sdp.primal))
Explanation: Solving it, we get:
End of explanation
time0 = time()
moments = get_moment_constraints(N, GHZ_state,measurements, GHZ_operators, 
                                 lambda_, 2)
extra = get_fullmeasurement([0 for _ in range(N)], measurements)
extramonomials = [extra]
moments[extra] = (1 - lambda_)
extra = get_fullmeasurement(flatten([1, [0 for _ in range(N-1)]]), measurements)
extramonomials.append(extra)
moments[extra] = (1/sqrt(2))*(1 - lambda_)
print("Constraints were generated in " + str(time()-time0) + " seconds.")
sdp = SdpRelaxation(flatten(measurements), parameters=[lambda_], verbose=1, parallel=True)
sdp.get_relaxation(2, objective=lambda_, substitutions=substitutions, 
                   momentsubstitutions=moments, extramonomials=extramonomials)
print("SDP relaxation was generated in " + str(time()-time0) + " seconds.")
sdp.solve(solver="mosek")
print("SDP was solved in " + str(sdp.solution_time) + " seconds.")
print("lambda_min is " + str(sdp.primal))
Explanation: Dual inequality
From the dual of the noise robustness SDP, one can extract a Bell inequality that detects the observed nonlocal correlations. The derived inequality will be in the form $\sum \alpha_{j_1 j_2 j_3 j_4}^{i_1 i_2 i_3 i_4} \langle M_{j_1}^{(i_1)}M_{j_2}^{(i_2)}M_{j_3}^{(i_3)}M_{j_4}^{(i_4)} \rangle + 1 \geq 0$.
We will present the code we used for the case of the $GHZ$ and $m = 2$ measurement choices. First, we solve the noise robustness SDP:
End of explanation
time0 = time()
sdp2 = SdpRelaxation(flatten(measurements), verbose=1, parallel=True)
sdp2.get_relaxation(2, substitutions=substitutions, extramonomials=extramonomials)
print("Second SDP relaxation was generated in " + str(time()-time0) + " seconds.")
Explanation: We only assign the values of up to the two-body correlators. This will ensure that we get the same expression for the Bell inequality as the one presented in the manuscript. 
Now, we proceed to generate a second relaxation, but we do not substitute the monomials corresponding to the observed correlations.
End of explanation
sdp2.status = sdp.status
sdp2.y_mat = [sdp.y_mat[1]]
Explanation: Then we assign to the second SDP the same dual variables as for the one solved before.
End of explanation
time0 = time()
bound = sdp2.extract_dual_value(1)
ineq = 0
for monomial, index in sdp2.monomial_index.items():
    ineq += round(2*sdp2.extract_dual_value(monomial),4)*monomial
print("Dual was generated in " + str(time()-time0) + " seconds")
Explanation: We extract the value of the dual variable for each monomial appearing in the moment matrix. For the function "extract_dual_value" to work properly we needed to generate a moment matrix in which the values of the correlators where not substituted yet.
End of explanation
print(ineq/bound)  
Explanation: Lastly, we normalize the inequality so to be in the form presented before (i.e. with the classical bound being equal to $1$ )
End of explanation
<END_TASK>