Compatibility_with_bipartite_models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computational appendix of [arXiv:1904.08943](https://arxiv.org/abs/1904.08943).\n",
    "---\n",
    "In this appendix we describe with explicit formulas the *scalar extension* method presented in [arXiv:1904.08943](https://arxiv.org/abs/1904.08943) to impose relaxations of (non-convex) factorization constraints in semidefinite programs, and show how this method can be used to study correlations in complex causal networks.\n",
    "\n",
    "authors: Alejandro Pozas-Kerstjens\n",
    "\n",
    "requires: [ncpol2sdpa](https://ncpol2sdpa.readthedocs.io/en/stable/), numpy, [picos](https://picos-api.gitlab.io/picos/) <= 1.1.3.post8, scipy, time\n",
    "\n",
    "last update: Apr, 2019\n",
    "\n",
    "Description of the method\n",
    "---\n",
    "Scalar extension is a method designed to deliver a hierarchy which, by means of semidefinite programming, allows to detect distributions that are not realizable in complex quantum networks. Essentially, the method consists in extending the generating set of monomials $\\mathcal{S}$ for a given moment matrix with additional monomials of the form $S_i\\langle S_j\\rangle,\\,S_i\\langle S_j\\rangle\\langle S_k\\rangle\\dots$ (throughout this notebook, the bracket and trace notation of expectation values will be used interchangeably). This is, we complement the generating set of monomials with another set, whose elements are products of operators multiplied by (possibly unknown) expectation values of other products of operators.\n",
    "\n",
    "This is particularly useful to test correlations in networks where some variables are independent of each other. An example of such is the tripartite-line scenario, depicted below\n",
    "\n",
    "<img src='tripartite_line.png' width='500' height='500'>\n",
    "\n",
    "A probability distribution $p(abc|xyz)$ is said to have a quantum realization in this scenario if there exists sets of projectors $\\{\\Pi_a^x\\}$, $\\{\\Pi_b^y\\}$ and $\\{\\Pi_c^z\\}$ (acting on $\\mathcal{H}_A$, $\\mathcal{H}_B=\\mathcal{H}_{B_1}\\otimes\\mathcal{H}_{B_2}$ and $\\mathcal{H}_C$, respectively), and a quantum state of the form $\\rho=\\rho_{AB_1}\\otimes\\rho_{B_2C}$ such that $p(abc|xyz)=\\text{Tr}(\\Pi_a^x\\otimes\\Pi_b^y\\otimes\\Pi_c^z\\rho)$. Moreover, from the NPA hierarchy it is known that if such a quantum realization exists, then for any set $\\mathcal{S}$ of products of projectors the matrix given by $\\Gamma_{ij}=\\text{Tr}(S_i^\\dagger S_j \\rho)$ is positive-semidefinite.\n",
    "\n",
    "However, a moment matrix generated from a distribution compatible with the tripartite-line scenario satisfies additional properties that matrices compatible with quantum realizations in other causal structures would not. Namely, if $S_i$ and $S_j$ only contain projectors of $A$ and $C$, the element in $\\Gamma_{ij}$ can be written as the product of two other elements in the matrix. This is, if $S_i^\\dagger S_j= S^A S^C$, where $S^A$ (resp. $S^C$) is a sequence of only projectors of $A$ (resp. $C$), then the moment matrix satisfies $\\Gamma_{ij}=\\Gamma_{a_1,a_2}\\Gamma_{c_1,c_2}$, where $S_{a_1}^\\dagger S_{a_2}=S^A$ and $S_{c_1}^\\dagger S_{c_2}=S^C$.\n",
    "\n",
    "Therefore, in order to know whether a distribution has a quantum realization in the tripartite-line scenario, we must be able to check whether moment matrices with unknown entries but that satisfy factorization constraints as that explained above can be made positive-semidefinite, this is, we must be able to solve problems of the form\n",
    "\n",
    "$$\\begin{aligned}\n",
    "& \\text{find}\n",
    "& & \\text{vars} \\\\\n",
    "& \\text{such that}\n",
    "& & \\Gamma \\geq 0.\n",
    "\\end{aligned}$$\n",
    "\n",
    "Failing to find a positive-semidefinite completion of one of these moment matrices will be a proof that the distribution does not have a quantum realization in the scenario. Nevertheless, factorization constraints are non-linear and furthermore non-convex, so they cannot be imposed in semidefinite programs. Scalar extension overcomes this issue by imposing linear constraints in larger moment matrices, that implement relaxations of the factorizations.\n",
    "\n",
    "Let us demonstrate the logic behind the method through an example, where we build a moment matrix step by step: assume we are given a binary probability distribution $p(abc|xyz)$ where $a,b,c,x,y,z\\in\\{0,1\\}$ and we want to test whether it has a quantum realization in the tripartite-line scenario. We will build a $4\\times 4$ moment matrix that, in principle, could allow us to answer the question in the negative: if such a matrix cannot be made positive-semidefinite, then the distribution does not have a quantum realization. For simplicity, we will assume that the entries in the matrix are real.\n",
    "\n",
    "Let us begin with an arbitrary, $4\\times 4$ symmetric matrix\n",
    "$$\n",
    "\\Gamma = \\begin{pmatrix}\n",
    "v_1 & v_2 & v_3 & v_4 \\\\\n",
    " & v_5 & v_6 & v_7 \\\\\n",
    " & & v_8 & v_9 \\\\\n",
    " & & & v_{10}\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "The first step proceeds as in the standard NPA hierarchy: we assume that a quantum realization exists and with it the set of projectors $\\{\\Pi_a^x\\}\\cup\\{\\Pi_b^y\\}\\cup\\{\\Pi_c^z\\}$. With it, we define a set of generating moments. In our case, let us choose\n",
    "\n",
    "$$\\mathcal{S}=\\{\\mathbb{1},\\Pi_{a=0}^{x=0},\\Pi_{a=0}^{x=0}\\Pi_{a=0}^{x=1},\\Pi_{c=0}^{z=0}\\Pi_{c=0}^{z=1}\\},$$\n",
    "\n",
    "and try to fill $\\Gamma$ with it according to $\\Gamma_{ij}=\\text{Tr}(S_i^\\dagger S_j \\rho)$. This allows us to reduce the number of variables in the problem. Since projectors are Hermitian and idempotent (so $\\Pi^\\dagger\\Pi=\\Pi^2=\\Pi\\,\\forall\\Pi$) we have that $v_2=\\text{Tr}(\\Pi_{a=0}^{x=0}\\rho)=\\text{Tr}({\\Pi_{a=0}^{x=0}}^\\dagger\\Pi_{a=0}^{x=0}\\rho)=v_5$. Proceeding analogously we also obtain $v_3=v_6$. We can also identify some of the elements of the matrix with elements of the probability distribution. Indeed, we have that $v_1=\\text{Tr}(\\rho)=1$ (which is true for any quantum state $\\rho$) and $v_2=\\text{Tr}(\\Pi_{a=0}^{x=0}\\rho)=p_a(0|0)$, where $p_A(a|x)$ is the marginalization of $p(abc|xyz)$ over the outcomes of parties $B$ and $C$. After all these substitutions, $\\Gamma$ has the form\n",
    "$$\n",
    "\\Gamma = \\begin{pmatrix}\n",
    "1 & p_A(0|0) & v_3 & v_4 \\\\\n",
    " & p_A(0|0) & v_3 & v_7 \\\\\n",
    " & & v_8 & v_9 \\\\\n",
    " & & & v_{10}\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "If the probability distribution was compatible with the tripartite-line scenario, expectation values of operator strings containing only projectors of $A$ and $C$ would factorize. This generates two additional constraints, namely\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "v_7 &= \\text{Tr}(\\Pi_{a=0}^{x=0}\\Pi_{c=0}^{z=0}\\Pi_{c=0}^{z=1}\\rho)=\\text{Tr}(\\Pi_{a=0}^{x=0}\\rho_{AB_1})\\text{Tr}(\\Pi_{c=0}^{z=0}\\Pi_{c=0}^{z=1}\\rho_{B_2C})=p_A(0|0)v_4,\\\\\n",
    "v_9 &= \\text{Tr}({\\Pi_{a=0}^{x=1}}^\\dagger{\\Pi_{a=0}^{x=0}}^\\dagger\\Pi_{c=0}^{z=0}\\Pi_{c=0}^{z=1}\\rho)=\\text{Tr}({\\Pi_{a=0}^{x=1}}^\\dagger{\\Pi_{a=0}^{x=0}}^\\dagger\\rho_{AB_1})\\text{Tr}(\\Pi_{c=0}^{z=0}\\Pi_{c=0}^{z=1}\\rho_{B_2C})=v_3 v_4.\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The first constraint can be readily imposed since it is a linear relation between $v_7$ and $v_4$ (recall that $p_A(0|0)$ is a known real number). Note that this is a nonlinear constraint, but it turns out that it can be imposed as a linear relation between unknown variables given a specific probability distribution. Therefore, this constraint can only be be imposed when aiming to certify if a correlation is compatible or not, and never for optimizing Bell inequalities over the set of compatible correlations.\n",
    "\n",
    "The second constraint involves products of two unknown variables, and therefore it can never be imposed in semidefinite programming. In order to impose a relaxation of such constraint, we perform the *scalar extension*. In particular, let us extend the moment matrix with just one new column (and the corresponding row), so the new moment matrix is\n",
    "\n",
    "$$\n",
    "\\tilde{\\Gamma} = \\begin{pmatrix}\n",
    "1 & p_A(0|0) & v_3 & v_4 & v_{11}\\\\\n",
    " & p_A(0|0) & v_3 & p_A(0|0)v_4 & v_{12} \\\\\n",
    " & & v_8 & v_9 & v_{13} \\\\\n",
    " & & & v_{10} & v_{14} \\\\\n",
    " & & & & v_{15}\n",
    "\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "We choose to associate such a new column with the operator $S_5=\\Pi_{a=0}^{x=1}\\left\\langle \\Pi_{c=0}^{z=1}\\Pi_{c=0}^{z=0}\\right\\rangle$, which is the projector $\\Pi_{a=0}^{x=1}$ multiplied by the scalar factor $\\left\\langle \\Pi_{c=0}^{z=1}\\Pi_{c=0}^{z=0}\\right\\rangle$, which cannot be computed from the probability distribution, and therefore will remain as a variable to be determined in our problem. Using as extended set of generating moments $\\tilde{\\mathcal{S}}=\\mathcal{S}\\cup\\{S_5\\}$, some new relations can be imposed in $\\tilde{\\Gamma}$, namely $v_{11}=p_a(0|1)v_4$ and $v_{14}=v_{15}$, obtaining\n",
    "\n",
    "$$\n",
    "\\tilde{\\Gamma} = \\begin{pmatrix}\n",
    "1 & p_A(0|0) & v_3 & v_4 & p_A(0|1)v_4\\\\\n",
    " & p_A(0|0) & v_3 & p_A(0|0)v_4 & v_{12} \\\\\n",
    " & & v_8 & v_9 & v_{13} \\\\\n",
    " & & & v_{10} & v_{14} \\\\\n",
    " & & & & v_{14}\n",
    "\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "Now, at last, the factorization condition translates to the linear constraint $v_9=v_{12}$. Therefore, in order to know whether the probability distribution $p(abc|xyz)$ has a quantum realization in the tripartite-line scenario, one can solve the problem\n",
    "\n",
    "\n",
    "$$\\begin{aligned}\n",
    "& \\text{find}\n",
    "& & v_3, v_4, v_8, v_9, v_{10}, v_{13}, v_{14} \\\\\n",
    "& \\text{such that}\n",
    "& & \\begin{pmatrix}\n",
    "1 & p_A(0|0) & v_3 & v_4 & p_A(0|1)v_4\\\\\n",
    " & p_A(0|0) & v_3 & p_A(0|0)v_4 & v_9 \\\\\n",
    " & & v_8 & v_9 & v_{13} \\\\\n",
    " & & & v_{10} & v_{14} \\\\\n",
    " & & & & v_{14}\n",
    "\\end{pmatrix} \\geq 0.\n",
    "\\end{aligned}$$\n",
    "\n",
    "This problem can be cast as a semidefinite program which, upon failing of finding a variable assignment that satisfied the constraints, would identify $p(abc|xyz)$ as not realizable. However, on the other hand success on finding such a variable assignment is not a proof of the distribution having a quantum realization, as in the standard NPA hierarchy.\n",
    "\n",
    "Application: discarding bi-quantum and bi-local models in the tripartite-line scenario\n",
    "---\n",
    "Now, we proceed to give an example of use of the method described above, by showing how it rightly identifies probability distributions that do not have realizations in the tripartite-line scenario in two different cases: when the sources distribute quantum states, and when they distribute classical variables. Explicitly, we show that the probability distribution $P^{22}$ in [Phys. Rev. A 85, 032119 (2012)](https://link.aps.org/doi/10.1103/PhysRevA.85.032119), when mixed with white noise giving rise to the noisy distribution $P_v=vP^{22}+(1-v)P_0$, fails to have a quantum realization in the tripartite-line scenario for $v>1/2$, and fails to have a bi-local hidden variable model for $v>1/4$.\n",
    "\n",
    "For doing so, we perform scalar extension to the moment matrices obtained in the standard application of the NPA hierarchy. However, instead of looking for a variable assignment that makes the moment matrix positive-semidefinite, due to better numerical stability we will solve the equivalent problem of optimizing the smallest eigenvalue of $\\tilde{\\Gamma}$. This is, we will solve problems of the form\n",
    "\n",
    "$$\\begin{aligned}\n",
    "& \\text{maximize}\n",
    "& & \\lambda \\\\\n",
    "& \\text{such that}\n",
    "& & \\tilde{\\Gamma} - \\lambda\\mathbb{1}\\geq 0.\n",
    "\\end{aligned}$$\n",
    "\n",
    "In this problem a positive optimal value of $\\lambda$ indicates that there exists a variable assignment for which $\\tilde{\\Gamma}$ is positive-semidefinite, and a negative optimal value identifies the distribution employed to construct $\\tilde{\\Gamma}$ as not compatible with the scenario under scrutiny. We will show that in the case of non-commuting operators the optimal value of $\\lambda$ is negative whenever $v > 1/2$, and whenever $v > 1/4$ when we additionally impose all operators to commute.\n",
    "\n",
    "To improve readability of this notebook and ease of use, we placed the supporting functions to a separate file; please download it 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](http://sdpa.sourceforge.net/) as an executable in the path or [Mosek](https://mosek.com/) with its Python interface installed; cvxopt as a solver is not recommended) together with the [Ncpol2sdpa](http://pypi.python.org/pypi/ncpol2sdpa) package.\n",
    "\n",
    "First, we import everything we will need:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scalar_extension_tools import fix_moments, get_factorization_constraints, get_moments_extracols\n",
    "from ncpol2sdpa import flatten, generate_variables, generate_operators,    \\\n",
    "                       projective_measurement_constraints, SdpRelaxation\n",
    "from numpy import eye\n",
    "from time import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the sake of simplicity, we define functions that generate sets of commuting and noncommuting variables (taken from [this example](https://github.com/FlavioBaccari/Hierarchy-for-nonlocality-detection)), a function that, given a probability distribution and a set of measurements, generates all computable quantities inside the moment matrix, and another function that will transform the feasibility problem into an optimization problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_commuting_measurements(party, label):\n",
    "    \"\"\"Generates the list of symbolic variables representing the measurements\n",
    "    for a given party. The variables are treated as commuting.\n",
    "\n",
    "    :param party: configuration indicating the configuration of number m\n",
    "                  of measurements and outcomes d for each measurement. It is a\n",
    "                  list with m integers, each of them representing the number of\n",
    "                  outcomes of the corresponding  measurement.\n",
    "    :type party: list of int\n",
    "    :param label: label to represent the given party\n",
    "    :type label: str\n",
    "\n",
    "    :returns: list of sympy.core.symbol.Symbol\n",
    "    \"\"\"\n",
    "    measurements = []\n",
    "    for i, p in enumerate(party):\n",
    "        measurements.append(generate_variables(label + '_%s' % i, p - 1,\n",
    "                                               hermitian=True))\n",
    "    return measurements\n",
    "\n",
    "def generate_noncommuting_measurements(party, label):\n",
    "    \"\"\"Generates the list of symbolic variables representing the measurements\n",
    "    for a given party. The variables are treated as noncommuting.\n",
    "\n",
    "    :param party: configuration indicating the configuration of number m\n",
    "                  of measurements and outcomes d for each measurement. It is a\n",
    "                  list with m integers, each of them representing the number of\n",
    "                  outcomes of the corresponding  measurement.\n",
    "    :type party: list of int\n",
    "    :param label: label to represent the given party\n",
    "    :type label: str\n",
    "\n",
    "    :returns: list of sympy.core.symbol.Symbol\n",
    "    \"\"\"\n",
    "    measurements = []\n",
    "    for i, p in enumerate(party):\n",
    "        measurements.append(generate_operators(label + '_%s' % i, p - 1,\n",
    "                                               hermitian=True))\n",
    "    return measurements\n",
    "\n",
    "def get_moment_constraints(probability, measurements):\n",
    "    \"\"\"Generates the list of moment equalities substitution for the SDP,\n",
    "    that is, the values of the correlators for a given probapility\n",
    "    distribution.\n",
    "    \n",
    "    :param probability: probability distribution.\n",
    "    :type probability: object\n",
    "    :param measurements: list of variables representing the local\n",
    "                         measurements in the given scenario.\n",
    "    :type measurements: list of lists of sympy.core.symbol.Symbol\n",
    "    \n",
    "    :returns: dict of sympy.core.add.Add\n",
    "    \"\"\"\n",
    "    measA, measB, measC = measurements\n",
    "    moments = {measA[0][0]: probability.p_a(0, 0),\n",
    "               measA[1][0]: probability.p_a(0, 1),\n",
    "               measB[0][0]: probability.p_b(0, 0),\n",
    "               measB[1][0]: probability.p_b(0, 1),\n",
    "               measC[0][0]: probability.p_c(0, 0),\n",
    "               measC[1][0]: probability.p_c(0, 1),\n",
    "               measA[0][0] * measB[0][0]: probability.p_ab(0, 0, 0, 0),\n",
    "               measA[0][0] * measB[1][0]: probability.p_ab(0, 0, 0, 1),\n",
    "               measA[1][0] * measB[0][0]: probability.p_ab(0, 0, 1, 0),\n",
    "               measA[0][0] * measB[1][0]: probability.p_ab(0, 0, 1, 1),\n",
    "               measA[0][0] * measC[0][0]: probability.p_ac(0, 0, 0, 0),\n",
    "               measA[0][0] * measC[1][0]: probability.p_ac(0, 0, 0, 1),\n",
    "               measA[1][0] * measC[0][0]: probability.p_ac(0, 0, 1, 0),\n",
    "               measA[1][0] * measC[1][0]: probability.p_ac(0, 0, 1, 1),\n",
    "               measB[0][0] * measC[0][0]: probability.p_bc(0, 0, 0, 0),\n",
    "               measB[0][0] * measC[1][0]: probability.p_bc(0, 0, 0, 1),\n",
    "               measB[1][0] * measC[0][0]: probability.p_bc(0, 0, 1, 0),\n",
    "               measB[1][0] * measC[1][0]: probability.p_bc(0, 0, 1, 1),\n",
    "               measB[0][0] * measC[0][0]: probability.p_bc(0, 0, 0, 0),\n",
    "               measA[0][0] * measB[0][0] * measC[0][0]: probability.p(0, 0, 0, 0, 0, 0),\n",
    "               measA[0][0] * measB[0][0] * measC[1][0]: probability.p(0, 0, 0, 0, 0, 1),\n",
    "               measA[0][0] * measB[1][0] * measC[0][0]: probability.p(0, 0, 0, 0, 1, 0),\n",
    "               measA[0][0] * measB[1][0] * measC[1][0]: probability.p(0, 0, 0, 0, 1, 1),\n",
    "               measA[1][0] * measB[0][0] * measC[0][0]: probability.p(0, 0, 0, 1, 0, 0),\n",
    "               measA[1][0] * measB[0][0] * measC[1][0]: probability.p(0, 0, 0, 1, 0, 1),\n",
    "               measA[1][0] * measB[1][0] * measC[0][0]: probability.p(0, 0, 0, 1, 1, 0),\n",
    "               measA[1][0] * measB[1][0] * measC[1][0]: probability.p(0, 0, 0, 1, 1, 1)\n",
    "               }\n",
    "    return moments\n",
    "\n",
    "def add_lambda(problem):\n",
    "    \"\"\"Transforms the feasibility problem\n",
    "    find      vars\n",
    "    such that M >= 0\n",
    "    into the optimization problem\n",
    "    maximize  lambda\n",
    "    such that M - lambda * I >= 0\n",
    "    \n",
    "    :param problem: initial feasibility problem\n",
    "    :type problem: picos.problem.Problem\n",
    "    \n",
    "    :returns: picos.problem.Problem\n",
    "    \"\"\"\n",
    "    picos_problem = problem.convert_to_picos()\n",
    "    X = picos_problem.get_variable('X')\n",
    "    λ = picos_problem.add_variable('λ')\n",
    "    picos_problem.set_objective('max', λ)\n",
    "    constraint = picos_problem.get_constraint(0)\n",
    "    constraint.Exp1 -=  λ * eye(X.size[0])\n",
    "    return picos_problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can focus on the problem at hand. First, we define the noisy version of the probability distribution $P^{22}$, which we will contrast against bi-local and bi-quantum hidden variable models, along with its marginals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class P22(object):\n",
    "    def __init__(self, visibility):\n",
    "        self.v = visibility\n",
    "\n",
    "    # Full distribution\n",
    "    def p(self, a, b, c, x, y, z):\n",
    "        return self.v * (1 + (-1) ** (a + b + c + x*y + y*z))/8 + (1 - self.v)/8\n",
    "\n",
    "    # Two-party marginals\n",
    "    def p_ab(self, a, b, x, y):\n",
    "        return sum([self.p(a, b, c, x, y, 0) for c in range(2)])\n",
    "    def p_bc(self, b, c, y, z):\n",
    "        return sum([self.p(a, b, c, 0, y, z) for a in range(2)])\n",
    "    def p_ac(self, a, c, x, z):\n",
    "        return sum([self.p(a, b, c, x, 0, z) for b in range(2)])\n",
    "\n",
    "    # Single-party marginals\n",
    "    def p_a(self, a, x):\n",
    "        return sum([self.p_ac(a, c, x, 0) for c in range(2)])\n",
    "    def p_b(self, b, y):\n",
    "        return sum([self.p_ab(a, b, 0, y) for a in range(2)])\n",
    "    def p_c(self, c, z):\n",
    "        return sum([self.p_ac(a, c, 0, z) for a in range(2)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we define the parameters relevant to our scenario: the number of measurements and outputs of each of the parties, and the level of the NPA hierarchy we will consider."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "level = 3      # Level of the NPA hierarchy we will use\n",
    "oA = [2, 2]    # Number of outputs of each of Alice's measurements\n",
    "oB = [2, 2]    # Number of outputs of each of Bob's measurements\n",
    "oC = [2, 2]    # Number of outputs of each of Charlie's measurements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From now on, we treat each of the cases separately.\n",
    "\n",
    "### Case 1: Non-existence of bi-quantum models for $v>1/2$\n",
    "---\n",
    "We begin defining the symbolic variables representing the operators, and the constraints they satisfy. In this case, operators representing different measurements of a same party do not commute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Alice's measurements\n",
    "measA = generate_noncommuting_measurements(oA, 'A')\n",
    "        \n",
    "# Bob's measurements\n",
    "measB = generate_noncommuting_measurements(oB, 'B')\n",
    "\n",
    "# Charlie's measurements\n",
    "measC = generate_noncommuting_measurements(oC, 'C')\n",
    "\n",
    "measurements = [measA, measB, measC]\n",
    "substitutions = projective_measurement_constraints(measurements)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we fix the visibility and create the noisy probability distribution $P_v$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v = 1/2 + 1e-5\n",
    "p_v = P22(v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After the distribution and the state and measurements are defined, one can begin calculating cells in the moment matrix. This is performed in three steps: first, we generate the computable moments in the $\\Gamma$ (the standard, NPA moment matrix) with ``get_moment_constraints()``"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "moments = get_moment_constraints(p_v, measurements)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Second, given the measurements for the parties that are causally-disconnected, we generate variables for the new columns that will be needed for the extension and the substitution rules for factorizing moments with ``get_factorization_constraints()``. Let us stress that if one wishes to obtain distribution-free constraints (and thus use the dual of the problem solution as a non-linear Bell inequality, for instance), then only the constraints that do not make use of explicit values of the probability distribution must be used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "causally_factorizing_moments, extra_columns, \\\n",
    "extra_column_names = get_factorization_constraints(parties=[measA, measC],\n",
    "                                                   moments=moments,\n",
    "                                                   substitutions=substitutions,\n",
    "                                                   level=level,\n",
    "                                                   return_column_names=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As described in the main text, we choose the new columns to be indexed by monomials of the form $\\left\\langle S_i \\right\\rangle\\mathbb{1}$, where $S_i$ is in the set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "extra_column_names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Third, for the new columns, we generate substitutions for the remaining moments that appear of which one of the factors is known (like $v_7=p_a(0|0) v_4$ in the example of the introduction) with ``get_moments_extracols()``. Again, these constraints should not be added in the case of looking for general certificates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "other_factorizing_moments = get_moments_extracols(moments, extra_columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before proceeding, we do a bit of post-processing in order to eliminate numerical-precision errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "moments = {**moments, **causally_factorizing_moments, **other_factorizing_moments}\n",
    "moments = fix_moments(moments, 1e-10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can generate the extended moment matrix by complementing the standard NPA moment matrix $\\Gamma$ at the level specified by `level` with the columns represented by `extra_columns`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time0 = time()\n",
    "sdp = SdpRelaxation(flatten(measurements))\n",
    "sdp.get_relaxation(level=level,                        # Initial level of the moment matrix\n",
    "                   extramonomials=extra_columns,       # Columns with which will be extended\n",
    "                   substitutions=substitutions,        # Idempotency and commutation rules\n",
    "                   momentsubstitutions=moments         # Variable substitutions\n",
    "                   )\n",
    "print(\"SDP relaxation was generated in \" + str(time()-time0) + \" seconds.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally we build solve the optimization problem\n",
    "$$\\begin{aligned}\n",
    "& \\text{maximize}\n",
    "& & \\lambda \\\\\n",
    "& \\text{such that}\n",
    "& & \\tilde{\\Gamma} - \\lambda\\mathbb{1}\\geq 0.\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "problem = add_lambda(sdp)\n",
    "problem.solve(solver='mosek', solve_via_dual=False, verbose=1, tol=1e-12)\n",
    "print('λ: ' + str(problem.obj_value()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "By varying the parameter $v$ it can be checked that for $v>1/2$ (up to numerical precision) the optimal value for the smallest eigenvalue of $\\tilde{\\Gamma}$ is negative, which is a certificate that the generated correlations cannot be explained by the sources sending quantum states and the parties performing measurements on their shares."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Case 2: Non-existence of bi-local models for $v>1/4$\n",
    "---\n",
    "As a last step, we additionally impose that variables representing different measurements of the same operator commute, and show that in this case no positive-semidefinite completion of the extended moment matrix can be found for $v>1/4$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Alice's measurements\n",
    "measA = generate_commuting_measurements(oA, 'A')\n",
    "subsA = projective_measurement_constraints(measA)\n",
    "\n",
    "# Bob's measurements\n",
    "measB = generate_commuting_measurements(oB, 'B')\n",
    "subsB = projective_measurement_constraints(measB)\n",
    "\n",
    "# Charlie's measurements\n",
    "measC = generate_commuting_measurements(oC, 'C')\n",
    "subsC = projective_measurement_constraints(measC)\n",
    "\n",
    "measurements = [measA, measB, measC]\n",
    "substitutions = {**subsA, **subsB, **subsC}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_v.v = 1/4 + 1e-5       # We set the new visibility here\n",
    "moments = get_moment_constraints(p_v, measurements)\n",
    "causally_factorizing_moments, extra_columns, \\\n",
    "extra_column_names = get_factorization_constraints(parties=[measA, measC],\n",
    "                                                   moments=moments,\n",
    "                                                   substitutions=substitutions,\n",
    "                                                   level=level,\n",
    "                                                   return_column_names=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the set of extra generating monomials is again $\\{\\left\\langle S_i \\right\\rangle\\mathbb{1}\\}_i$, where $S_i$ is in"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "extra_column_names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us continue with the calculations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time0 = time()\n",
    "other_factorizing_moments = get_moments_extracols(moments, extra_columns)\n",
    "\n",
    "moments = {**moments, **causally_factorizing_moments, **other_factorizing_moments}\n",
    "moments = fix_moments(moments, 1e-10)\n",
    "print(\"Constraints were generated in \" + str(time()-time0) + \" seconds.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time0 = time()\n",
    "sdp = SdpRelaxation(flatten(measurements))\n",
    "sdp.get_relaxation(level=level, extramonomials=extra_columns, substitutions=substitutions, momentsubstitutions=moments)\n",
    "print(\"SDP relaxation was generated in \" + str(time()-time0) + \" seconds.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "probl = add_lambda(sdp)\n",
    "probl.solve(solver='mosek', solve_via_dual=False, verbose=1, tol=1e-12)\n",
    "print('λ: ' + str(probl.obj_value()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indeed, the smallest eigenvalue of $\\tilde{\\Gamma}$, is negative for $v > 1/4$, which certifies that such distributions cannot be explained with bi-local hidden variable models."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}