!pip install qiskit[visualization]

# Imports for Qiskit
from qiskit import QuantumCircuit, execute, Aer, IBMQ
from qiskit.compiler import transpile, assemble
from qiskit.tools.jupyter import *
from qiskit.visualization import *
from qiskit import *
from qiskit import IBMQ, Aer
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister
from qiskit.visualization import plot_histogram
import numpy as np

# Various imports 
import numpy as np

from copy import deepcopy
from matplotlib import pyplot as plt

# IBMQ.save_account('Put your token')
# provider = IBMQ.load_account()
# IBMQ.get_provider(hub='ibm-q', group='open', project = 'main')

# Create the various registers needed
clock = QuantumRegister(2, name='clock')
input = QuantumRegister(1, name='b')
ancilla = QuantumRegister(1, name='ancilla')
measurement = ClassicalRegister(2, name='c')

# Create an empty circuit with the specified registers
circuit = QuantumCircuit(ancilla, clock, input, measurement)

circuit.barrier()
circuit.draw(output='mpl')

def qft_dagger(circ, q, n):      
    circ.h(clock[1]);
    for j in reversed(range(n)):
      for k in reversed(range(j+1,n)):
        circ.cu1(-np.pi/float(2**(k-j)), q[k], q[j]);
    circ.h(clock[0]);
    circ.swap(clock[0], clock[1]);

def qft(circ, q, n):
    circ.swap(clock[0], clock[1]);
    circ.h(clock[0]);
    for j in reversed(range(n)):
      for k in reversed(range(j+1,n)):
        circ.cu1(np.pi/float(2**(k-j)), q[k], q[j]);
    circ.h(clock[1]);

def qpe(circ, clock, target):
    circuit.barrier()

    # e^{i*A*t}
    circuit.cu(np.pi/2, -np.pi/2, np.pi/2, 3*np.pi/4, clock[0], input, label='U');
    
    # e^{i*A*t*2}
    circuit.cu(np.pi, np.pi, 0, 0, clock[1], input, label='U2');

    circuit.barrier();
    
    # Perform an inverse QFT on the register holding the eigenvalues
    qft_dagger(circuit, clock, 2)
    
def inv_qpe(circ, clock, target):
    
    # Perform a QFT on the register holding the eigenvalues
    qft(circuit, clock, 2)

    circuit.barrier()

    # e^{i*A*t*2}
    circuit.cu(np.pi, np.pi, 0, 0, clock[1], input, label='U2');

    #circuit.barrier();

    # e^{i*A*t}
    circuit.cu(np.pi/2, np.pi/2, -np.pi/2, -3*np.pi/4, clock[0], input, label='U');

    circuit.barrier()



def hhl(circ, ancilla, clock, input, measurement):
    
    qpe(circ, clock, input)

    circuit.barrier()
    
    # This section is to test and implement C = 1
    circuit.cry(np.pi, clock[0], ancilla)
    circuit.cry(np.pi/3, clock[1], ancilla)

    circuit.barrier()
    
    circuit.measure(ancilla, measurement[0])
    circuit.barrier()
    inv_qpe(circ, clock, input)

# State preparation. (various initial values)
# (restart from second cell if changed)

intial_state = [0,1]
# intial_state = [1,0]
# intial_state = [1/np.sqrt(2),1/np.sqrt(2)]
# intial_state = [np.sqrt(0.9),np.sqrt(0.1)]

circuit.initialize(intial_state, 3)

circuit.barrier()

# Perform a Hadamard Transform
circuit.h(clock)

hhl(circuit, ancilla, clock, input, measurement)

# Perform a Hadamard Transform
circuit.h(clock)

circuit.barrier()


circuit.measure(input, measurement[1])

circuit.draw('mpl',scale=1)
#print(circuit)

# Execute the circuit using the simulator
simulator = qiskit.BasicAer.get_backend('qasm_simulator')
job = execute(circuit, backend=simulator, shots=1000)

#Get the result of the execution
result = job.result()

# Get the counts, the frequency of each answer
counts = result.get_counts(circuit)

# Display the results
plot_histogram(counts)

bcknd = Aer.get_backend('statevector_simulator')

job_sim = execute(circuit, bcknd)
result = job_sim.result()

o_state_result = result.get_statevector(circuit, decimals=3)
print(o_state_result)

provider.backends()

# Choose the backend on which to run the circuit
backend = provider.get_backend('ibmq_santiago')

from qiskit.tools.monitor import job_monitor

# Execute the job
job_exp = execute(circuit, backend=backend, shots=8192)

# Monitor the job to know where we are in the queue
job_monitor(job_exp, interval = 2)

# Get the results from the computation
results = job_exp.result()

# Get the statistics
answer = results.get_counts(circuit)

# Plot the results
plot_histogram(answer)



# Auto generated circuit (almost) matching the form from the Wong paper (using cu gates)

from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister

beta = 0
cycle_time = 0.5

qr = QuantumRegister(4)
cr = ClassicalRegister(4)
qc = QuantumCircuit(qr, cr, name="main")

intial_state = [0,1]
qc.initialize(intial_state, qr[3])

qc.barrier(qr[0], qr[1], qr[2], qr[3])

qc.h(qr[1])
qc.h(qr[2])

qc.barrier(qr[0], qr[1], qr[2], qr[3])

# e^{i*A*t}
qc.cu(np.pi/2, -np.pi/2, np.pi/2, 3*np.pi/4, qr[1], qr[3], label='U');

# e^{i*A*t*2}
qc.cu(np.pi, np.pi, 0, 0, qr[2], qr[3], label='U2');

qc.barrier(qr[0], qr[1], qr[2], qr[3])
qc.h(qr[2])
qc.cu1(-3.14159/2, qr[1], qr[2])
qc.h(qr[1])
qc.swap(qr[2], qr[1])
qc.barrier(qr[0], qr[1], qr[2], qr[3])
qc.cry(3.14159, qr[1], qr[0])
qc.cry(3.14159/3, qr[2], qr[0])

qc.barrier(qr[0], qr[1], qr[2], qr[3])
qc.swap(qr[2], qr[1])
qc.h(qr[1])
qc.cu1(3.14159/2, qr[1], qr[2])
qc.h(qr[2])
qc.barrier(qr[0], qr[1], qr[2], qr[3])

# e^{i*A*t*2}
qc.cu(np.pi, np.pi, 0, 0, qr[2], qr[3], label='U2');

# e^{i*A*t}
qc.cu(np.pi/2, np.pi/2, -np.pi/2, -3*np.pi/4, qr[1], qr[3], label='U');

qc.barrier(qr[0], qr[1], qr[2], qr[3])
qc.h(qr[1])
qc.h(qr[2])

qc.barrier(qr[0], qr[1], qr[2], qr[3])
qc.measure(qr[0], cr[0])
qc.measure(qr[3], cr[3])

from qiskit import execute, Aer
backend = Aer.get_backend("qasm_simulator")  # Use Aer qasm_simulator

job = execute(qc, backend, shots=1000)
result = job.result()
counts = result.get_counts(qc)
print("Total counts are:", counts)

# Draw the circuit
print(qc)
#qc.draw('mpl',scale=1)

# Plot a histogram
from qiskit.visualization import plot_histogram
plot_histogram(counts)


qc.draw('mpl',scale=1)




# Same example as above but using the U3 and U1 gates instead of U1
# Initialize with RY instead of "initialize()"

from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister

beta = 1

qr = QuantumRegister(4)
cr = ClassicalRegister(4)
qc = QuantumCircuit(qr, cr, name="main")

# intial_state = [0,1]
# qc.initialize(intial_state, qr[3])
qc.ry(3.14159*beta, qr[3])

qc.barrier(qr[0], qr[1], qr[2], qr[3])

qc.h(qr[1])
qc.h(qr[2])
qc.barrier(qr[0], qr[1], qr[2], qr[3])

# e^{i*A*t}
# qc.cu(np.pi/2, -np.pi/2, np.pi/2, 3*np.pi/4, qr[1], qr[3], label='U');

# The CU gate is equivalent to a CU1 on the control bit followed by a CU3
qc.u1(3*3.14159/4, qr[1])
qc.cu3(3.14159/2, -3.14159/2, 3.14159/2, qr[1], qr[3])

# e^{i*A*t*2}
# qc.cu(np.pi, np.pi, 0, 0, qr[2], qr[3], label='U2');
qc.cu3(3.14159, 3.14159, 0, qr[2], qr[3])

qc.barrier(qr[0], qr[1], qr[2], qr[3])
qc.h(qr[2])
qc.cu1(-3.14159/2, qr[1], qr[2])
qc.h(qr[1])
qc.swap(qr[2], qr[1])
qc.barrier(qr[0], qr[1], qr[2], qr[3])
qc.cry(3.14159, qr[1], qr[0])
qc.cry(3.14159/3, qr[2], qr[0])
qc.barrier(qr[0], qr[1], qr[2], qr[3])
qc.swap(qr[2], qr[1])
qc.h(qr[1])
qc.cu1(3.14159/2, qr[1], qr[2])
qc.h(qr[2])
qc.barrier(qr[0], qr[1], qr[2], qr[3])

# e^{i*A*t*2}
# qc.cu(np.pi, np.pi, 0, 0, qr[2], qr[3], label='U2');
qc.cu3(3.14159, 3.14159, 0, qr[2], qr[3])

# e^{i*A*t}
# qc.cu(np.pi/2, np.pi/2, -np.pi/2, -3*np.pi/4, qr[1], qr[3], label='U');

# The CU gate is equivalent to a CU1 on the control bit follwed by a CU3
qc.u1(-3*3.14159/4, qr[1])
qc.cu3(3.14159/2, 3.14159/2, -3.14159/2, qr[1], qr[3])

qc.barrier(qr[0], qr[1], qr[2], qr[3])

qc.h(qr[1])
qc.h(qr[2])

qc.barrier(qr[0], qr[1], qr[2], qr[3])

qc.measure(qr[0], cr[0])
qc.measure(qr[3], cr[3])
# qc.measure(qr, cr)

from qiskit import execute, Aer
backend = Aer.get_backend("qasm_simulator")  # Use Aer qasm_simulator

job = execute(qc, backend, shots=1000)
result = job.result()
counts = result.get_counts(qc)
print("Total counts are:", counts)

# Draw the circuit
print(qc)

# Plot a histogram
from qiskit.visualization import plot_histogram
plot_histogram(counts)


qc.draw('mpl',scale=1)


