In [1]:
% load_ext autoreload
% autoreload 2
% matplotlib inline

from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
import scipy as sp
import seaborn as sns
np.set_printoptions(precision=5, linewidth=120)
from tqdm import *
from drift_qec.A import *


/Users/yan/.miniconda/lib/python2.7/site-packages/matplotlib/__init__.py:872: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.
  warnings.warn(self.msg_depr % (key, alt_key))

In [18]:
r = np.linspace(-1.0, 1.0, 21)
b1, b2, b3, b4, b5, b6 = np.meshgrid(r,r,r,r,r,r)
b1 = np.ravel(b1)
b2 = np.ravel(b2)
b3 = np.ravel(b3)
b4 = np.ravel(b4)
b5 = np.ravel(b5)
b6 = np.ravel(b6)
check = (b1 ** 2 + b2 ** 2 + b3 ** 2 + b4 ** 2 + b5 ** 2 + b6 ** 2) <= 1
b1 = b1[check].copy()
b2 = b2[check].copy()
b3 = b3[check].copy()
b4 = b4[check].copy()
b5 = b5[check].copy()
b6 = b6[check].copy()

In [23]:
def p1L1(b1, b2, b3, b4, b5, b6, d):
    return (b1*np.cos(d) - b4*np.sin(d)) ** 2 + (b2*np.sin(d)) ** 2
def p2L1(b1, b2, b3, b4, b5, b6, d):
    return (b1*np.sin(d) + b4*np.cos(d)) ** 2 + (b2*np.cos(d)) ** 2
def p3L1(b1, b2, b3, b4, b5, b6, d):
    return b3 ** 2 + b5 ** 2 + b6 ** 2

def p1L2(b1, b2, b3, b4, b5, b6, d):
    return (b1*np.cos(d) - b5*np.sin(d)) ** 2 + (b3 ** 2 + b6 ** 2)*np.sin(d) ** 2
def p2L2(b1, b2, b3, b4, b5, b6, d):
    return b2 ** 2 + b4 ** 2
def p3L2(b1, b2, b3, b4, b5, b6, d):
    return (b1*np.sin(d) + b5*np.cos(d)) ** 2 + (b3 ** 2 + b6 ** 2)*np.cos(d) ** 2

def p1L3(b1, b2, b3, b4, b5, b6, d):
    return b1 ** 2
def p2L3(b1, b2, b3, b4, b5, b6, d):
    return (b2*np.cos(d) - b6*np.sin(d)) ** 2 + (b4*np.cos(d) - b5*np.sin(d)) ** 2 + (b3 * np.sin(d)) ** 2
def p3L3(b1, b2, b3, b4, b5, b6, d):
    return (b2*np.sin(d) + b6*np.cos(d)) ** 2 + (b4*np.sin(d) + b5*np.cos(d)) ** 2 + (b3 * np.cos(d)) ** 2

def update_prior(x, b1, b2, b3, b4, b5, b6, d, priorB):
    update = x[0] * np.log(p1L1(b1, b2, b3, b4, b5, b6, d)) \
           + x[1] * np.log(p2L1(b1, b2, b3, b4, b5, b6, d)) \
           + x[2] * np.log(p3L1(b1, b2, b3, b4, b5, b6, d)) \
           + x[3] * np.log(p1L2(b1, b2, b3, b4, b5, b6, d)) \
           + x[4] * np.log(p2L2(b1, b2, b3, b4, b5, b6, d)) \
           + x[5] * np.log(p3L2(b1, b2, b3, b4, b5, b6, d)) \
           + x[6] * np.log(p1L3(b1, b2, b3, b4, b5, b6, d)) \
           + x[7] * np.log(p2L3(b1, b2, b3, b4, b5, b6, d)) \
           + x[8] * np.log(p3L3(b1, b2, b3, b4, b5, b6, d)) 
    new_prior = update + priorB
    return new_prior

def reconstruct_M(b1, b2, b3, b4, b5, b6, idx):
    c1, c2, c3, c4, c5, c6 = b1[idx], b2[idx], b3[idx], b4[idx], b5[idx], b6[idx]
    B = np.array([
            [c1, c4, c5],
            [ 0, c2, c6],
            [ 0,  0, c3]
        ])
    M = np.dot(B.T, B)
    return M

In [52]:
D = 0.1
channel = Channel(kx=0.7, ky=0.2, kz=0.1, Q=FIXEDQ,
                  n=1e3, d1=D, d2=D, d3=D)
priorB = np.ones(b1.shape) * (-np.log(float(np.prod(b1.shape))))

In [56]:
data = channel.sample_data()
print data
priorB = update_prior(data, b1, b2, b3, b4, b5, b6, D, priorB)

idx = np.argmax(new_prior)
channel.M = reconstruct_M(b1, b2, b3, b4, b5, b6, idx)
channel.M


[ 0.5    0.343  0.157  0.508  0.337  0.155  0.539  0.316  0.145]
Out[56]:
array([[ 0.49,  0.35,  0.35],
       [ 0.35,  0.25,  0.25],
       [ 0.35,  0.25,  0.26]])

In [51]:
np.max(priorB)


Out[51]:
-33.847920125492116

In [50]:
np.mean(priorB > -100)


Out[50]:
0.88809640981743543

In [47]:
np.mean(priorB == -np.inf)


Out[47]:
0.11116946115972685

In [26]:
channel.C


Out[26]:
array([[ 0.13174,  0.11094, -0.02851],
       [ 0.11094,  0.55535, -0.21218],
       [-0.02851, -0.21218,  0.31285]])

In [ ]: