from qiskit import *
from qiskit.tools.visualization import plot_histogram
import numpy as np

qr = QuantumRegister(1)
cr = ClassicalRegister(1)

error = {}
for n in range(1,11):

    # Create a blank circuit
    qc = QuantumCircuit(qr,cr)
    
    # Implement an approximate Hadamard
    theta = np.pi/np.sqrt(2) # here we correctly choose theta=pi/sqrt(2)
    for j in range(n):
        qc.rx(theta/n,qr[0])
        qc.rz(theta/n,qr[0])
      
    # We need to measure how good the above approximation is. Here's a simple way to do this.
    # Step 1: Use a real hadamard to cancel the above approximation.
    # For a good approximatuon, the qubit will return to state 0. For a bad one, it will end up as some superposition.
    qc.h(qr[0])
    
    # Step 2: Run the circuit, and see how many times we get the outcome 1.
    # Since it should return 0 with certainty, the fraction of 1s is a measure of the error.
    qc.measure(qr,cr)
    shots = 20000
    job = execute(qc, Aer.get_backend('qasm_simulator'),shots=shots)
    try:
        error[n] = (job.result().get_counts()['1']/shots)
    except:
        pass
        
plot_histogram(error)

# The linear nature of error^(-1/2) shows that the error has a quadratic decay.
inverse_square_of_error = {}
for n in error:
    inverse_square_of_error[n] = (error[n])**(-1/2)
plot_histogram(inverse_square_of_error)

qr = QuantumRegister(1)
cr = ClassicalRegister(1)

error = {}
for n in range(1,11):

    # Create a blank circuit
    qc = QuantumCircuit(qr,cr)
    
    # Implement an approximate Hadamard
    theta = np.pi/np.sqrt(2) # here we correctly use theta=pi/sqrt(2)
    for j in range(n):
        qc.rz(theta/(2*n),qr[0])
        qc.rx(theta/n,qr[0])
        qc.rz(theta/(2*n),qr[0])
      
    # We need to measure how good the above approximation is. Here's a simple way to do this.
    # Step 1: Use a real hadamard to cancel the above approximation.
    # For a good approximatuon, the qubit will return to state 0. For a bad one, it will end up as some superposition.
    qc.h(qr[0])
    
    # Step 2: Run the circuit, and see how many times we get the outcome 1.
    # Since it should return 0 with certainty, the fraction of 1s is a measure of the error.
    qc.measure(qr,cr)
    shots = 100000
    job = execute(qc, Aer.get_backend('qasm_simulator'),shots=shots)
    try:
        error[n] = (job.result().get_counts()['1']/shots)
    except:
        pass
        
plot_histogram(error)

# The linear nature of error^(-1/3) shows that the error has a cubic decay.
# Note: this needs loads of shots to get a good result.
inverse_cube_of_error = {}
for n in error:
    error[n]
    inverse_cube_of_error[n] = (error[n])**(-1/3)
plot_histogram(inverse_cube_of_error)

