From 124edad8e34389f7aa319db4fcb6474523879dc4 Mon Sep 17 00:00:00 2001
From: Daniel Tsvetkov <danieltcv@gmail.com>
Date: Fri, 27 Mar 2020 18:59:07 +0100
Subject: [PATCH] own simulator krisi2

---
 09_quantum_monthy.py |  3 ++
 lib.py               | 73 ++++++++++++++++++++++++++++++++++++--------
 2 files changed, 64 insertions(+), 12 deletions(-)

diff --git a/09_quantum_monthy.py b/09_quantum_monthy.py
index 533a24f..5278b7d 100644
--- a/09_quantum_monthy.py
+++ b/09_quantum_monthy.py
@@ -55,9 +55,11 @@ def perform_exp(case):
 
 if __name__ == "__main__":
     correct = 0
+    results = defaultdict(int)
     for i in range(100):
         case = random.choice([CASE_IDENTICAL, CASE_ORTHOGONAL])
         result = perform_exp(case=case)
+        results[result] += 1
         if result == '11':
             guess = CASE_ORTHOGONAL
         else:
@@ -65,3 +67,4 @@ if __name__ == "__main__":
         if guess == case:
             correct += 1
     print("Correct: {}".format(correct))
+    # print("Results: {}".format(results))
diff --git a/lib.py b/lib.py
index 9a32782..77b7567 100644
--- a/lib.py
+++ b/lib.py
@@ -172,18 +172,9 @@ class State(Vector):
     def _is_normalized(self):
         return np.isclose(np.sum(np.abs(self.m ** 2)), 1.0)
 
-    @staticmethod
-    def _normalize_angles(theta, phi):
-        # theta is between [0 and pi]
-        theta = theta.real % np.pi if theta.real != np.pi else theta.real
-        # phi is between (0 and 2*pi]
-        phi = phi.real % (2 * np.pi)
-        return theta, phi
-
     @classmethod
     def from_bloch_angles(cls, theta, phi):
         """Creates a state from angles in Bloch sphere"""
-        theta, phi = cls._normalize_angles(theta, phi)
         m0 = np.cos(theta / 2)
         m1 = np.sin(theta / 2) * np.power(np.e, (1j * phi))
         m = m0 * _0 + m1 * _1
@@ -845,7 +836,7 @@ class MeasurementOperator(HermitianOperator):
         # <ψ|
         state_ct = state.conjugate_transpose()
         #  This is: <ψ|   .   M†_m M_m . |ψ>
-        return state_ct.m.dot(self.m.dot(state.m)).item()
+        return state_ct.m.dot(self.m.dot(state.m)).item().real
 
 
 m0 = MeasurementOperator.create_from_basis(_0)
@@ -1304,7 +1295,63 @@ def test_light():
     experiment(3, random_lights, [ver_filter, diag_filter, hor_filter])
 
 
-def test_krisi_measurement():
+def test_krisi_measurement_2():
+    from qiskit import (
+        QuantumCircuit,
+        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>")))
@@ -1325,4 +1372,6 @@ if __name__ == "__main__":
     test()
     test_quantum_processor()
     test_light()
-    test_krisi_measurement()
+    test_krisi_measurement_2()
+    # TODO:
+    test_krisi_measurement_3()