Finish example on quantum steering

[peterwittek]

The hierachy for quantum steering is supported, but a working example exampel is missing. The core idea is simple, it follows the NPA hierarchy:

from ncpol2sdpa import SteeringHierarchy, convert_to_picos, flatten,\
                       generate_measurements, \
                       projective_measurement_constraints,

nA, nB, mA, mB, dA, dB, dC = 3, 3, 2, 2, 2, 2, 2  # System configuration
level = 1  # Level of relaxation

A = generate_measurements([mA for _ in range(nA)], 'A')
B = generate_measurements([mB for _ in range(nB)], 'B')
substitutions = projective_measurement_constraints(A, B)
steering = SteeringHierarchy(flatten([A, B]), matrix_var_dim=dC, verbose=2)
steering.get_relaxation(level, substitutions=substitutions,
                        extramonomials=[Ai*Bj for Ai in flatten(A)
                                        for Bj in flatten(B)])
P = convert_to_picos(steering)

From here, we have to specify the assemblage and the rest of the constraints in PICOS.

[peterwittek]

The assemblage would be something like this:

from numpy import array, sqrt, hstack, kron
from qutip import tensor, sigmax, sigmay, sigmaz, qeye, Qobj, ptrace

w = 1  # Noise level in W state
IdC = qeye(dC)

def get_assemblage():
    # W state
    psi = Qobj(array([0, 1, 1, 0, 1, 0, 0, 0])/sqrt(3))
    rho = Qobj(w*(psi*psi.dag()).data+(1-w)*qeye(8).data/8,
               dims=[[2, 2, 2], [2, 2, 2]])
    # Measurement choices
    A = [[] for _ in range(nA)]
    A[0].append((qeye(2)+sigmax())/2)
    A[0].append(qeye(2)-A[0][-1])
    A[1].append((qeye(2)+sigmay())/2)
    A[1].append(qeye(2)-A[1][-1])
    A[2].append((qeye(2)+sigmaz())/2)
    A[2].append(qeye(2)-A[2][-1])
    B = [[] for _ in range(nB)]
    B[0].append((qeye(2)+sigmax())/2)
    B[0].append(qeye(2)-B[0][-1])
    B[1].append((qeye(2)+sigmay())/2)
    B[1].append(qeye(2)-B[1][-1])
    B[2].append((qeye(2)+sigmaz())/2)
    B[2].append(qeye(2)-B[2][-1])

    # Assemblage
    assemblage = [[[[[] for _ in range(nA)] for _ in range(nB)]
                   for _ in range(mA)] for _ in range(mB)]
    for x in range(nA):
        for y in range(nB):
            for a in range(mA):
                for b in range(mB):
                    assemblage[a][b][x][y] = \
                      ptrace(tensor(A[x][a], B[y][b], IdC)*rho, 2)
    return assemblage