#import Aer here, before calling qiskit_ionq_provider
from qiskit import Aer

from qiskit_ionq import IonQProvider 

#Call provider and set token value
provider = IonQProvider(token='my token')

provider.backends()

from qiskit import *
#import qiskit.aqua as aqua
#from qiskit.quantum_info import Pauli
#from qiskit.aqua.operators.primitive_ops import PauliOp
from qiskit.circuit.library import PhaseEstimation
from qiskit import QuantumCircuit

from matplotlib import pyplot as plt
import numpy as np


from qiskit.chemistry.drivers import PySCFDriver, UnitsType
from qiskit.chemistry.drivers import Molecule

def buildControlledT(p, m):

    # initialize the circuit
    qc = QuantumCircuit(2, 1)

    # Hardmard on ancilla, now in |+>
    qc.h(0)

    # initialize to |1>
    qc.x(1)

    # applying T gate to qubit 1
    for i in range(2**m):
        qc.cp(np.pi/4, 0, 1)

    # phase correction
    qc.rz(p, 0)

    # inverse QFT (in other words, just measuring in the x-basis)
    qc.h(0)

    qc.measure([0],[0])

    return qc

def IPEA(k, backend_string):

    # get backend
    if backend_string == 'qpu':
        backend = provider.get_backend('ionq_qpu')
    elif backend_string == 'qasm':
        backend = Aer.get_backend('qasm_simulator')

    # bits
    bits = []

    # phase correction
    phase = 0.0
    
    # loop over iterations 
    for i in range(k-1, -1, -1):
       
        # construct the circuit
        qc = buildControlledT(phase, i)
         
        # run the circuit
        job = execute(qc, backend)
        
        if backend_string == 'qpu':
            from qiskit.providers.jobstatus import JobStatus
            import time
            # Check if job is done
            while job.status() is not JobStatus.DONE:
                print("Job status is", job.status() )
                time.sleep(60)

            # grab a coffee! This can take up to a few minutes.
            # once we break out of that while loop, we know our job is finished
            print("Job status is", job.status() )
            print(job.get_counts()) # these counts are the “true” counts from the actual QPU Run
        
        # get result
        result = job.result()

        # get current bit
        this_bit = int(max(result.get_counts(), key=result.get_counts().get))
        print(result.get_counts())

        bits.append(this_bit)

        # update phase correction
        phase /= 2
        phase -= (2 * np.pi * this_bit / 4.0)

    return bits

def eig_from_bits(bits):
    
    eig = 0.

    m = len(bits)

    # loop over all bits
    for k in range(len(bits)):
        eig += bits[k] / (2**(m-k))

    #eig *= 2*np.pi

    return eig

# perform IPEA
backend = 'qasm'
bits = IPEA(5, backend)
print(bits)

# re-construct energy
eig = eig_from_bits(bits)

print(eig)

#perform IPEA with different values of n
n_values = []
eig_values = []
for i in range(1, 8):
    n_values.append(i)
    # perform IPEA
    backend = 'qasm'
    bits = IPEA(i, backend)

    # re-construct energy
    eig = eig_from_bits(bits)
    eig_values.append(eig)

n_values, eig_values = np.array(n_values), np.array(eig_values)
plt.plot(n_values, eig_values)
plt.xlabel('n (bits)', fontsize=15)
plt.ylabel(r'$\phi$', fontsize=15)
plt.title(r'$\phi$ vs. n', fontsize=15)

# perform IPEA
backend = 'qpu'
bits = IPEA(5, backend)
print(bits)

# re-construct energy
eig = eig_from_bits(bits)

print(eig)

#perform IPEA with different values of n
n_values = []
eig_values = []
for i in range(1, 8):
    n_values.append(i)
    # perform IPEA
    backend = 'qpu'
    bits = IPEA(i, backend)

    # re-construct energy
    eig = eig_from_bits(bits)
    eig_values.append(eig)

n_values, eig_values = np.array(n_values), np.array(eig_values)
plt.plot(n_values, eig_values)
plt.xlabel('n (bits)', fontsize=15)
plt.ylabel(r'$\phi$', fontsize=15)
plt.title(r'$\phi$ vs. n', fontsize=15)