In [1]:

    

%load_ext autoreload
%autoreload 2
%matplotlib inline

from matplotlib import pyplot as plt
from cvxopt import matrix, solvers
solvers.options['show_progress'] = False
import numpy as np
import pandas as pd
import scipy as sp
import seaborn as sns
np.set_printoptions(precision=3, linewidth=120)
from copy import copy
from tqdm import *
from drift_qec.Q 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 [9]:

    

matrix?


    

In [3]:

    

D = 0.01
channel = Channel(kx=0.7, ky=0.2, kz=0.1,
                  Q=np.linalg.qr(np.random.randn(3,3))[0],
                  n=1e3, d1=D, d2=D, d3=D)
data = channel.sample_data()


    

In [55]:

    

L = SENSOR(D,D,D)
Ld = L[:, [0, 3, 5]]
m11, m22, m33 = list(np.dot(np.dot(np.linalg.inv(np.dot(Ld.T, Ld)), Ld.T), data))
m12 = (0.5*(m11-m22)*np.cos(2.0*D) - 0.5*(data[0]-data[1])) / np.sin(2.0*D)
m13 = (0.5*(m11-m33)*np.cos(2.0*D) - 0.5*(data[3]-data[5])) / np.sin(2.0*D)
m23 = (0.5*(m22-m33)*np.cos(2.0*D) - 0.5*(data[7]-data[8])) / np.sin(2.0*D)
b = np.array([[m12 ** 2 / np.sqrt(m22 * m33)],
              [m13 ** 2 / np.sqrt(m11 * m33)],
              [m23 ** 2 / np.sqrt(m11 * m22)]])
x = b / np.linalg.norm(b)
x = np.array([np.sqrt(x[0] * np.sqrt(m22 * m33)),
              np.sqrt(x[1] * np.sqrt(m11 * m33)),
              np.sqrt(x[2] * np.sqrt(m11 * m22))])
x = x.T[0]

x = np.sign([m12, m13, m23]).T * x
Mhat = np.array([
        [ m11, x[0], x[1]],
        [x[0],  m22, x[2]],
        [x[1], x[2], m33]
    ])


    

In [59]:

    

Mhat


    

    
Out[59]:





array([[ 0.467,  0.489, -0.349],
       [ 0.489,  0.18 ,  0.168],
       [-0.349,  0.168,  0.353]])

In [60]:

    

np.linalg.eig(Mhat)


    

    
Out[60]:





(array([-0.358,  0.905,  0.453]), array([[-0.588,  0.808,  0.044],
        [ 0.674,  0.459,  0.579],
        [-0.448, -0.371,  0.814]]))

In [61]:

    

channel.C


    

    
Out[61]:





array([[ 0.462,  0.105, -0.229],
       [ 0.105,  0.188, -0.144],
       [-0.229, -0.144,  0.351]])

In [25]:

    

b


    

    
Out[25]:





matrix([[ 2.412],
        [-1.068],
        [ 0.72 ]])

In [24]:

    

A = matrix(np.eye(3))
b = matrix([[m12 / np.sqrt(m22 * m33)],
            [m13 / np.sqrt(m11 * m33)],
            [m23 / np.sqrt(m11 * m22)]])
G = matrix([[-1, 0, 0, 0, 1, 0, 0],
            [ 0,-1, 0, 0, 0, 1, 0],
            [ 0, 0,-1, 0, 0, 0, 1]])
h = matrix( [ 0, 0, 0, 1, 0, 0, 0])
x = solvers.coneqp(A.T*A, -A.T*b, G, h)['x']
print(x)


    

    




---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-24-28e206c44ff5> in <module>()
      7             [ 0, 0,-1, 0, 0, 0, 1]])
      8 h = matrix( [ 0, 0, 0, 1, 0, 0, 0])
----> 9 x = solvers.coneqp(A.T*A, -A.T*b, G, h)['x']
     10 print(x)

/Users/yan/.miniconda/lib/python2.7/site-packages/cvxopt/coneprog.pyc in coneqp(P, q, G, h, dims, A, b, initvals, kktsolver, xnewcopy, xdot, xaxpy, xscal, ynewcopy, ydot, yaxpy, yscal, **kwargs)
   1820     if (not matrixP or (not matrixG and G is not None) or 
   1821         (not matrixA and A is not None)) and not customkkt:
-> 1822         raise ValueError("use of function valued P, G, A requires a "\
   1823             "user-provided kktsolver")
   1824     customx = (xnewcopy != None or xdot != None or xaxpy != None or 

ValueError: use of function valued P, G, A requires a user-provided kktsolver

In [ ]:

    

x = Variable(3)
Q = matrix([[m11,0.0,0.0], [0.0,m22,0.0], [0.0,0.0,m33]])
c = m11*m22*m33
objective = Minimize( (x[0]-m12) ** 2 + (x[1] - m13) ** 2 + (x[2] - m23) ** 2 )
constraints = [ quad_form(x,Q) <= c ]

p = Problem(objective, constraints)
primal_result = p.solve()
x = np.array(x.value).T[0]


    

In [ ]:

    

objective.value


    

In [ ]:

    

Mhat = np.array([
        [m11,  x[0], x[1]],
        [x[0],  m22, x[2]],
        [x[1], x[2],  m33],
    ])


    

In [ ]:

    

Mhat


    

In [ ]:

    

np.trace(Mhat)


    

In [ ]:

    

channel.C


    

In [ ]:

    

np.linalg.eig(Mhat)[0]


    

In [ ]: