using Convex

function get_visibility(K)
    noise = real([trace(K[i])*eye(2)/2 for i=1:size(K, 1)])
    P = [[ComplexVariable(2, 2) for i=1:2] for j=1:6]
    q = Variable(6, Positive())
    t = Variable(1, Positive())
    constraints = [P[i][j] in :SDP for i=1:6 for j=1:2]
    constraints += sum(q)==1
    constraints += t<=1
    constraints += [P[i][1]+P[i][2] == q[i]*eye(2) for i=1:6]
    constraints += t*K[1] + (1-t)*noise[1] == P[1][1] + P[2][1] + P[3][1]
    constraints += t*K[2] + (1-t)*noise[2] == P[1][2] + P[4][1] + P[5][1]
    constraints += t*K[3] + (1-t)*noise[3] == P[2][2] + P[4][2] + P[6][1]
    constraints += t*K[4] + (1-t)*noise[4] == P[3][2] + P[5][2] + P[6][2]
    p = maximize(t, constraints)
    solve!(p)
    return p.optval
end

function dp(v)
    eye(2) + v[1]*[0 1; 1 0] + v[2]*[0 -im; im 0] + v[3]*[1 0; 0 -1]
end
b = [ 1  1  1; 
     -1 -1  1; 
     -1  1 -1;  
      1 -1 -1]/sqrt(3)
M = [dp(b[i, :]) for i=1:size(b,1)]/4;
get_visibility(M)