Meaningful conjugation code

abp/clifford.py

@@ -15,9 +15,9 @@ decompositions = ("xxxx", "xx", "zzxx", "zz", "zxx", "z", "zzz", "xxz",
                   "xzx", "xzxxx", "xzzzx", "xxxzx", "xzz", "zzx", "xxx", "x",
                   "zzzx", "xxzx", "zx", "zxxx", "xxxz", "xzzz", "xz", "xzxx")
 
-def conjugate(vop, transform):
+def conjugate(operator, unitary):
     """ Returns transform * vop * transform^dagger and a phase in {+1, -1} """
-    return measurement_table[vop, transform]
+    return measurement_table[operator, unitary]
 
 def use_old_cz():
     """ Use the CZ table from A&B's code """
@@ -101,22 +101,32 @@ def get_state_table(unitaries):
             np.dot(kp, state).T)
     return state_table
 
+def get_measurement_entry(operator, unitary):
+    """ 
+    Any Clifford group unitary will map an operator A in {I, X, Y, Z} 
+    to an operator B in +-{I, X, Y, Z}. This finds that mapping.
+    """
+    matrices = ({"x": qi.msqx, "z": qi.sqz}[c] for c in decompositions[unitary])
+    unitary = reduce(np.dot, matrices, np.eye(2, dtype=complex))
+    operator = qi.operators[operator]
+    new_operator = reduce(np.dot, (unitary, operator, qi.hermitian_conjugate(unitary)))
+
+    for i, o in enumerate(qi.operators):
+        if np.allclose(o, new_operator):
+            return i, 1
+        elif np.allclose(o, -new_operator):
+            return i, -1
+
+    raise IndexError
+
 def get_measurement_table():
     """ 
     Compute a table of transform * operation * transform^dagger 
     This is pretty unintelligible right now, we should probably compute the phase from unitaries instead
     """
     measurement_table = np.zeros((4, 24, 2), dtype=complex)
-    for vop, transform in it.product(range(4), range(24)):
-        assert vop in set(xrange(4))
-        op = times_table[transform, vop]
-        op = times_table[op, conjugation_table[transform]]
-        is_id_or_vop = (transform % 4 == 0) or (transform % 4 == vop)
-        is_non_pauli = (transform >= 4) and (transform <= 15)
-        phase = ((-1, 1), (1, -1))[is_id_or_vop][is_non_pauli]
-        if vop == 0: 
-            phase = 1
-        measurement_table[vop, transform] = [op, phase]
+    for operator, unitary in it.product(range(4), range(24)):
+        measurement_table[operator, unitary] = [operator, unitary]
     return measurement_table
 
 

abp/qi.py

@@ -36,6 +36,8 @@ cz = np.array(np.eye(4), dtype=complex)
 cz[3,3]=-1
 
 # States
+zero = np.array([[1],[0]], dtype=complex)
+one = np.array([[0],[1]], dtype=complex)
 plus = np.array([[1],[1]], dtype=complex) / np.sqrt(2)
 bond = cz.dot(np.kron(plus, plus))
 nobond = np.kron(plus, plus)
@@ -46,6 +48,7 @@ names = "identity", "px", "py", "pz", "hadamard", "phase", "sqz", "msqz", "sqy",
 by_name = dict(zip(names, common_us))
 
 paulis = px, py, pz
+operators = id, px, py, pz
 
 
 def normalize_global_phase(m):

scripts/wtf.py

@@ -0,0 +1,46 @@
+from abp import clifford, qi
+import numpy as np
+import itertools as it
+
+def conj_hack(operator, unitary):
+    """ TODO: make this meaningful / legible """
+    assert operator in set(xrange(4))
+    op = clifford.times_table[unitary, operator]
+    op = clifford.times_table[op, clifford.conjugation_table[unitary]]
+    is_id_or_operator = (unitary % 4 == 0) or (unitary % 4 == operator)
+    is_non_pauli = (unitary >= 4) and (unitary <= 15)
+    phase = ((-1, 1), (1, -1))[is_id_or_operator][is_non_pauli]
+    if operator == 0: 
+        phase = 1
+    return op, phase
+
+def conj(operator, unitary):
+    """ Better """
+    matrices = ({"x": qi.msqx, "z": qi.sqz}[c] for c in clifford.decompositions[unitary])
+    unitary = reduce(np.dot, matrices, np.eye(2, dtype=complex))
+    operator = qi.operators[operator]
+    new_operator = reduce(np.dot, (unitary, operator, qi.hermitian_conjugate(unitary)))
+
+    for i, o in enumerate(qi.operators):
+        if np.allclose(o, new_operator):
+            return i, 1
+        elif np.allclose(o, -new_operator):
+            return i, -1
+
+    raise IndexError
+
+
+for operator, unitary in it.product(range(4), range(24)):
+    assert conj(operator, unitary) == conj_hack(operator, unitary)
+
+    #new = np.dot(u, np.dot(o, qi.hermitian_conjugate(u)))
+    #which = clifford.find_clifford(new, clifford.unitaries[:4])
+    #assert which in xrange(4)
+    #whichm = [qi.id, qi.px, qi.py, qi.pz][which]
+    #if np.allclose(new, whichm):
+        #return which, 1
+    #elif np.allclose(new, -whichm):
+        #return which, -1
+    #else:
+        #raise IndexError
+

tests/test_conjugation.py

@@ -4,7 +4,6 @@ import itertools
 
 #//! replaces op by trans * op * trans^dagger and returns a phase,
 #/*! either +1 or -1 (as RightPhase(0) or RightPhase(2)) */
-#RightPhase conjugate (const LocCliffOp trans);
 
 def test_conjugation():
     """ Test that clifford.conugate() agrees with graphsim.LocCliffOp.conjugate """