### calculate the probability for theta = 1 degree

import math 
import cmath 

def probs_unitary(ang):
  ang = math.radians(ang) # 1 degree
  y = complex(0, ang) #i1
  y0 = 0.5*(1 + cmath.exp(y))
  #print ('probability for getting 0: ', abs(y0)**2)
  y1 = 0.5*(1 - cmath.exp(y))
  #print ('probability for getting 1: ', abs(y1)**2)
  return abs(y0)**2, abs(y1)**2

prob0, prob1 = probs_unitary(1)
print (r'probability for getting 0 and 1 for $\theta=1$: ', prob0, prob1)
prob0_10, prob1_10 = probs_unitary(10)
print (r'probability for getting 0 and 1 for $\theta=10$: ', prob0_10, prob1_10)

from google.colab import  drive
drive.mount('/content/drive/')

!pip install qiskit

import qiskit as q


import matplotlib
import matplotlib.pyplot as plt
import numpy as np
# from mpl_toolkits.mplot3d import proj3d

print (q.__version__)

print(matplotlib.__version__)

style = {'backgroundcolor': 'lavender', 'dpi':200, 'subfontsize':10}

PEcircuit = q.QuantumCircuit(3, 2)
PEcircuit.x(2) # this is the eigen state of S gate; in figures above this $\psi$
PEcircuit.draw()

# start by applying the H gates
for qubit in range(2):
    PEcircuit.h(qubit)
PEcircuit.draw()

repetitions = 1
for counting_qubit in range(2):
    for i in range(repetitions):
        PEcircuit.cp(np.pi/2, counting_qubit, 2); # This is CU # QuantumCircuit.cp(theta, control_qubit, target_qubit, label=None, ctrl_state=None)
        # importantly to apply controlled phase rotation we need to know the rotation angle beforehand
    repetitions *= 2 # 2^{n-1} factor
PEcircuit.draw()

!pip install pylatexenc

style = {'backgroundcolor': 'lavender', 'dpi':200, 'subfontsize':10}
PEcircuit.draw(output='mpl', scale=0.8, style=style, filename='/content/drive/My Drive/Colab Notebooks/Quantum_Compute/FourierTransform/2bits_QPE_QiskitEx.png')

def qft_dagger(qc, n):
    """n-qubit inverse qft
       first n qubits will be transformed
       check the qft circuit implementation first 
    """
    for qubit in range(n//2):
        qc.swap(qubit, n-qubit-1)
    for j in range(n):
        for m in range(j):
            qc.cp(-math.pi/float(2**(j-m)), m, j)
        qc.h(j)

PEcircuit.barrier()
# Apply inverse QFT
qft_dagger(PEcircuit, 2)
# Measure
PEcircuit.barrier()
for n in range(2):
    PEcircuit.measure(n,n)



# PEcircuit.draw(output='mpl', style=style, filename='/content/drive/My Drive/Colab Notebooks/Quantum_Compute/FourierTransform/2bits_QPE_QiskitEx_Full.png')    

aer_sim = q.Aer.get_backend('aer_simulator')
shots = 2048
t_PEcircuit = q.transpile(PEcircuit, aer_sim)
qobj = q.assemble(t_PEcircuit, shots=shots)
results = aer_sim.run(qobj).result()
answer = results.get_counts()

q.visualization.plot_histogram(answer)

PEcircuit_1 = q.QuantumCircuit(4, 3)
PEcircuit_1.x(3) # this is the eigen state of S gate; in figures above this $\psi$

for qubit in range(3):
    PEcircuit_1.h(qubit)

repetitions = 1
for counting_qubit in range(3):
    for i in range(repetitions):
        PEcircuit_1.cp(np.pi/2, counting_qubit, 3) # This is CU # QuantumCircuit.cp(theta, control_qubit, target_qubit, label=None, ctrl_state=None)
        # importantly to apply controlled phase rotation we need to know the rotation angle beforehand
    repetitions *= 2 # 2^{n-1} factor
PEcircuit_1.draw()

PEcircuit_1.barrier()
# Apply inverse QFT
qft_dagger(PEcircuit_1, 3)
# Measure
PEcircuit_1.barrier()
for n in range(3):
    PEcircuit_1.measure(n, n)


PEcircuit_1.draw(output='mpl', style=style, filename='/content/drive/My Drive/Colab Notebooks/Quantum_Compute/FourierTransform/3bits_QPE_QiskitEx.png')    

aer_sim = q.Aer.get_backend('aer_simulator')
shots = 2048
t_PEcircuit_1 = q.transpile(PEcircuit_1, aer_sim)
qobj = q.assemble(t_PEcircuit_1, shots=shots)
results_1 = aer_sim.run(qobj).result()
answer_1 = results_1.get_counts()

q.visualization.plot_histogram(answer_1)

PEcircuit_1.num_qubits-1

def qpe_circ(nq, angle):
  
  """
  generic qpe circuit
  nq: number of qubits
  angle: theta as defined in description
  """

  qc = q.QuantumCircuit(nq+1, nq) # num qubits + eigen state, for measuring we need num qubits = classical bits for measure. 

  qc.x(nq) # apply the not gate to create state 1, i.e. eigen state at the last qubit
  for qubit in range(nq):
    qc.h(qubit)

  repetitions = 1
  for counting_qubit in range(qc.num_qubits - 1): #qc.num_qubits - 1 = nq 
    for i in range(repetitions):
        qc.cp(np.pi *2 * angle, counting_qubit, qc.num_qubits-1); # This is CU # QuantumCircuit.cp(theta, control_qubit, target_qubit, label=None, ctrl_state=None)
        # importantly to apply controlled phase rotation we need to know the rotation angle beforehand
    repetitions *= 2 # 2^{n-1} factor
  qc.barrier()
  qft_dagger(qc, qc.num_qubits-1)
  qc.barrier()
  for n1 in range(qc.num_qubits -1):
    qc.measure(n1, n1)
  return qc

circuit_1 = qpe_circ(3, 1/5)

circuit_1.draw(output='mpl', style=style, filename='/content/drive/My Drive/Colab Notebooks/Quantum_Compute/FourierTransform/3bits_QPE_QiskitEx_NonInt_Theta.png')

aer_sim = q.Aer.get_backend('aer_simulator')
shots = 2048
t_circuit_1 = q.transpile(circuit_1, aer_sim)
qobj = q.assemble(t_circuit_1, shots=shots)
results_c1 = aer_sim.run(qobj).result()
answer_c1 = results_c1.get_counts()

q.visualization.plot_histogram(answer_c1)

circuit_2 = qpe_circ(5, 1/5)

circuit_2.draw(output='mpl', style=style, filename='/content/drive/My Drive/Colab Notebooks/Quantum_Compute/FourierTransform/5bits_QPE_QiskitEx_NonInt_Theta.png')

aer_sim = q.Aer.get_backend('aer_simulator')
shots = 2048
t_circuit_2 = q.transpile(circuit_2, aer_sim)
qobj = q.assemble(t_circuit_2, shots=shots)
results_c2 = aer_sim.run(qobj).result()
answer_c2 = results_c2.get_counts()

fig, ax = plt.subplots(figsize=(14, 6))
q.visualization.plot_histogram(answer_c2, figsize=(14, 6), color='violet', title=r'QPE Measurement Result with 5 Qubits: $\theta = 1/5$', ax=ax)
fig.savefig('/content/drive/My Drive/Colab Notebooks/Quantum_Compute/FourierTransform/5bits_QPE_QiskitEx_NonInt_Theta_Histogram.png', dpi=250)

for k, v  in answer_c2.items():
  print(k, v)

