import numpy as np
from qiskit import (
    QuantumCircuit,
    execute,
    Aer)
from qiskit.visualization import plot_histogram, plot_bloch_vector
from matplotlib import pyplot as plt
from collections import defaultdict

# Use Aer's qasm_simulator
simulator = Aer.get_backend('qasm_simulator')


def create_bloch_vector(theta, phi):
    x = np.sin(theta) * np.cos(phi)
    y = np.sin(theta) * np.sin(phi)
    z = np.cos(theta)
    return [x, y, z]


def perform_exp(iterations, case, draw=False):
    simulator = Aer.get_backend('qasm_simulator')
    all_counts = defaultdict(int)

    for i in range(iterations):
        qc = QuantumCircuit(2, 2)

        # uniform distribution on the sphere
        t = round(np.random.uniform(0, 1), 10)
        theta = np.arccos(1 - 2 * t)
        phi = round(np.random.uniform(0, 2 * np.pi), 10)

        CASE_ANGLES = {
            'identical': (theta, phi),
            'orthogonal': (theta + np.pi, phi)
        }

        qc.r(theta, phi, 0)
        theta1, phi1 = CASE_ANGLES[case]
        qc.r(theta1, phi1, 1)

        # print(theta, phi, phi1)

        if draw:
            q0 = create_bloch_vector(theta, phi)
            plot_bloch_vector(q0)
            plt.show()
            q1 = create_bloch_vector(theta1, phi1)
            plot_bloch_vector(q1)
            plt.show()

        # Measure in the Bell's basis
        qc.cx(0, 1)
        qc.h(0)
        qc.measure([0, 1], [0, 1])

        # print(qc)

        # execute the job
        job = execute(qc, simulator, shots=1)
        result = job.result()
        counts = result.get_counts(qc)

        for k, v in counts.items():
            all_counts[k] += v

    # gather all counts
    print("\nTotal counts for experiment '{}' are:\n".format(case))
    for k, v in sorted(all_counts.items(), key=lambda x: x[0]):
        print("{}: {}".format(k, v))


if __name__ == "__main__":
    # draw produces a Bloch sphere but does not work in this notebook

    perform_exp(iterations=100, case='identical', draw=False)  # same state
    perform_exp(iterations=100, case='orthogonal', draw=False)  # different orthogonal states