attempt at 3qubits game

10_krisi_3qubits_game.py

@@ -0,0 +1,136 @@
+import random
+from collections import defaultdict
+
+import numpy as np
+
+from lib import State, b_phi_p, b_phi_m, b_psi_p, b_psi_m, s, generate_bins
+
+
+def test_krisi_measurement_2():
+    CASE_IDENTICAL = "identical"
+    CASE_ORTHOGONAL = "orthogonal"
+
+    def perform_exp(case):
+        # produce angles with uniform distribution on the sphere
+        t = round(np.random.uniform(0, 1), 10)
+        theta0 = np.arccos(1 - 2 * t)
+        phi0 = round(np.random.uniform(0, 2 * np.pi), 10)
+
+        # rotate the 0th qubit in (theta, phi)
+        q0 = State.from_bloch_angles(theta0, phi0)
+
+        # rotate the 1st qubit depending on the case we are exploring...
+        if case == CASE_IDENTICAL:
+            # ... for identical we are rotating the 1st qubit the same angles as
+            # the 0th
+            theta1, phi1 = theta0, phi0
+        else:
+            # for orthogonal we rotate the 1st qubit in 90 degrees
+            # orthogonal to the first
+            theta1, phi1 = theta0 + np.pi, phi0
+
+        q1 = State.from_bloch_angles(theta1, phi1)
+
+        # Measure in the Bell's basis
+        st = State(q0 * q1)
+        meas = st.measure(basis=[b_phi_p, b_phi_m, b_psi_p, b_psi_m])
+        return meas
+
+    correct = 0
+    cases = defaultdict(int)
+    results = defaultdict(int)
+    for i in range(1000):
+        case = random.choice([CASE_IDENTICAL, CASE_ORTHOGONAL])
+        cases[case] += 1
+        result = perform_exp(case=case)
+        results[result] += 1
+        if result == '11':
+            guess = CASE_ORTHOGONAL
+            assert guess == case
+        else:
+            guess = CASE_IDENTICAL
+        if guess == case:
+            correct += 1
+    print("Correct: {}".format(correct))
+    # print("Results: {}".format(results))
+    # print("Cases: {}".format(cases))
+
+
+def test_krisi_measurement_3():
+    CASE_ID_ID = "id_id"
+    CASE_ID_ORT = "id_ort"
+    CASE_ORT_ID = "ort_id"
+    CASE_ORT_ORT = "ort_ort"
+
+    beta = 0
+    b_0 = State(1 / np.sqrt(2) * (s("|100>") - s("|010>")))
+    b_1 = State(1 / np.sqrt(2) * (s("|011>") - s("|101>")))
+    basis = [
+        s("|000>"),
+        State(np.cos(beta) * b_0 + np.sin(beta) * s("|001>")),
+        State(-np.sin(beta) * b_0 + np.cos(beta) * s("|001>")),
+        b_0,
+        b_1,
+        State(np.cos(beta) * b_1 + np.sin(beta) * s("|001>")),
+        State(-np.sin(beta) * b_1 + np.cos(beta) * s("|001>")),
+        s("|111>"),
+    ]
+
+    def perform_exp(case):
+        # produce angles with uniform distribution on the sphere
+        t = round(np.random.uniform(0, 1), 10)
+        theta0 = np.arccos(1 - 2 * t)
+        phi0 = round(np.random.uniform(0, 2 * np.pi), 10)
+
+        # rotate the 0th qubit in (theta, phi)
+        q0 = State.from_bloch_angles(theta0, phi0)
+
+        # rotate the 1st qubit depending on the case we are exploring...
+        if case == CASE_ID_ID:
+            theta1, phi1 = theta0, phi0
+            theta2, phi2 = theta0, phi0
+        elif case == CASE_ID_ORT:
+            theta1, phi1 = theta0, phi0
+            theta2, phi2 = theta0 + np.pi, phi0
+        elif case == CASE_ORT_ID:
+            theta1, phi1 = theta0 + np.pi, phi0
+            theta2, phi2 = theta0, phi0
+        else:
+            # CASE_ORT_ORT
+            theta1, phi1 = theta0 + np.pi, phi0
+            theta2, phi2 = theta0 + np.pi, phi0
+
+        q1 = State.from_bloch_angles(theta1, phi1)
+        q2 = State.from_bloch_angles(theta2, phi2)
+
+        # Measure in the arbitrary basis
+        st = State(q0 * q1 * q2)
+        meas = st.measure(basis=basis)
+
+        return meas
+
+    def krisi_3_format_results(results):
+        all_pos = generate_bins(8)
+        print("Results:")
+        for case, rv in results.items():
+            print("{}:".format(case))
+            print("   raw  : {}".format(sorted(rv.items())))
+            case_total = sum(rv.values())
+            print("   total: {}".format(case_total))
+            for pos in all_pos:
+                value = rv.get(pos, 0)
+                percent = value / case_total
+                print("   {}  : {:3d} ({:.2f}%)".format(pos, value, percent))
+
+    results = defaultdict(lambda: defaultdict(int))
+    for i in range(1000):
+        case = random.choice([CASE_ID_ID, CASE_ID_ORT, CASE_ORT_ID, CASE_ORT_ORT])
+        result = perform_exp(case=case)
+        results[case][result] += 1
+
+    krisi_3_format_results(results)
+
+
+if __name__ == "__main__":
+    # test_krisi_measurement_2()
+    test_krisi_measurement_3()

lib.py

@@ -2,7 +2,6 @@ import random
 from collections import defaultdict
 from functools import reduce
 from numbers import Complex
-from pprint import pprint
 from typing import Union
 
 import numpy as np
@@ -273,7 +272,7 @@ class State(Vector):
     def get_fmt_of_element(self):
         return "{:0" + str(int(np.ceil(np.log2(len(self))))) + "b}"
 
-    def measure(self, basis=None):
+    def measure(self, basis=None, allow_empty=False):
         """
         m_ops: a set of MeasurementOperators.
                e.g. for computational basis, the m_ops is [|0><0|, |1><1|]
@@ -310,10 +309,18 @@ class State(Vector):
             prob = self.get_prob_from_measurement_op(m_op)
             weights += [prob]
 
+        empty_choices, empty_weights = [], []
+        if allow_empty == True:
             empty_prob = 1.0 - sum(weights)
+            empty_weights = [empty_prob]
+            empty_choices = [REPR_EMPTY_SET]
+        else:
+            # normalize the weights
+            weights = list(np.array(weights) / sum(weights))
+
         format_str = self.get_fmt_of_element()
-        choices = [REPR_EMPTY_SET] + [format_str.format(i) for i in range(len(weights))]
+        choices = empty_choices + [format_str.format(i) for i in range(len(weights))]
-        weights = [empty_prob] + weights
+        weights = empty_weights + weights
         self.measurement_result = random.choices(choices, weights)[0]
 
         return self.measurement_result
@@ -1278,7 +1285,7 @@ def test_light():
             for filter in filters:
                 st = filter.on(st)
             st = State(st)
-            total_light[st.measure()] += 1
+            total_light[st.measure(allow_empty=True)] += 1
         print("{}. {}:\n    {}".format(id, human_state, dict(total_light)))
 
     its = 100
@@ -1295,83 +1302,12 @@ def test_light():
     experiment(3, random_lights, [ver_filter, diag_filter, hor_filter])
 
 
-def test_krisi_measurement_2():
+def generate_bins(count):
-    from qiskit import (
+    format_str = "{:0" + str(int(np.ceil(np.log2(count)))) + "b}"
-        QuantumCircuit,
+    return [format_str.format(i) for i in range(count)]
-        execute,
-        Aer)
-    simulator = Aer.get_backend('qasm_simulator')
-
-    CASE_IDENTICAL = "identical"
-    CASE_ORTHOGONAL = "orthogonal"
-
-    def perform_exp(case):
-        # produce angles with uniform distribution on the sphere
-        t = round(np.random.uniform(0, 1), 10)
-        theta0 = np.arccos(1 - 2 * t)
-        phi0 = round(np.random.uniform(0, 2 * np.pi), 10)
-
-        # rotate the 0th qubit in (theta, phi)
-        q0 = State.from_bloch_angles(theta0, phi0)
-
-        # rotate the 1st qubit depending on the case we are exploring...
-        if case == CASE_IDENTICAL:
-            # ... for identical we are rotating the 1st qubit the same angles as
-            # the 0th
-            theta1, phi1 = theta0, phi0
-        else:
-            # for orthogonal we rotate the 1st qubit in 90 degrees
-            # orthogonal to the first
-            theta1, phi1 = theta0 + np.pi, phi0
-
-        q1 = State.from_bloch_angles(theta1, phi1)
-
-        # Measure in the Bell's basis
-        st = State(q0 * q1)
-        meas = st.measure(basis=[b_phi_p, b_phi_m, b_psi_p, b_psi_m])
-        return meas
-
-    correct = 0
-    cases = defaultdict(int)
-    results = defaultdict(int)
-    for i in range(1000):
-        case = random.choice([CASE_IDENTICAL, CASE_ORTHOGONAL])
-        cases[case] += 1
-        result = perform_exp(case=case)
-        results[result] += 1
-        if result == '11':
-            guess = CASE_ORTHOGONAL
-            assert guess == case
-        else:
-            guess = CASE_IDENTICAL
-        if guess == case:
-            correct += 1
-    print("Correct: {}".format(correct))
-    # print("Results: {}".format(results))
-    # print("Cases: {}".format(cases))
-
-
-def test_krisi_measurement_3():
-    beta = 2.6
-    b_0 = State(1 / np.sqrt(2) * (s("|100>") - s("|010>")))
-    b_1 = State(1 / np.sqrt(2) * (s("|011>") - s("|101>")))
-    meas = s("|000>").measure(basis=[
-        s("|000>"),
-        State(np.cos(beta) * b_0 + np.sin(beta) * s("|001>")),
-        State(-np.sin(beta) * b_0 + np.cos(beta) * s("|001>")),
-        b_0,
-        b_1,
-        State(np.cos(beta) * b_1 + np.sin(beta) * s("|001>")),
-        State(-np.sin(beta) * b_1 + np.cos(beta) * s("|001>")),
-        s("|111>"),
-    ])
-    print(meas)
 
 
 if __name__ == "__main__":
     test()
     test_quantum_processor()
     test_light()
-    test_krisi_measurement_2()
-    # TODO:
-    test_krisi_measurement_3()