import qiskit
from qiskit import QuantumCircuit
from qiskit.visualization import plot_histogram

import numpy as np

from qiskit import execute, Aer

simulator = Aer.get_backend('qasm_simulator')

import matplotlib.pyplot as plt

import warnings
warnings.filterwarnings('ignore')

from qiskit.quantum_info import DensityMatrix, Statevector, Operator

def getDensityMatrix(circuit):
    return DensityMatrix(circuit).data

state_0 = np.array([1, 0])
state_1 = np.array([0, 1])

from functools import reduce

Dag = lambda matrix: matrix.conj().T
Kron = lambda *matrices: reduce(np.kron, matrices)

def pm(matrix):
    for row in range(len(matrix)):
        for col in range (len(matrix[row])):
            print("{:.3f}".format(matrix[row][col]), end = " ")
        print()

def getMeasurements(n):
    psi_0 = np.array([1.0, 0.0])
    psi_1 = np.array([0.0, 1.0])
    I = np.eye(2)

    M_0 = psi_0.reshape([2, 1]) @ psi_0.reshape([1, 2]).conj()
    M_1 = psi_1.reshape([2, 1]) @ psi_1.reshape([1, 2]).conj()
    
    M = [M_0, M_1]
    
    measurements = []
    for i in range(2 ** n):
        binnum = bin(i)[2:].rjust(n, '0')

        temp = []
        indices = map(lambda x: int(x), list(binnum))
        for index in indices:
            temp.append(M[index])
        
        measurements.append( Kron(*temp) )
    
    return measurements

def measure(M, state):
    return Dag(state) @ Dag(M) @ M @ state

# https://qiskit.org/textbook/ch-algorithms/grover.html

def initCircuit(n):
    circuit = QuantumCircuit(n, n)
    
    for i in range(n):
        circuit.h(i)
    circuit.barrier()
    
    return circuit

inputCircuit_2q = initCircuit(2)
inputCircuit_2q.draw(output='mpl')

def createOracle_3():
    circuit = QuantumCircuit(2, 2)

    # Oracle for find 3

    # U_f 
    circuit.cz(0, 1)
    circuit.barrier()
    
    return circuit

oracleCircuit_3 = createOracle_3()
oracleCircuit_3.draw(output='mpl')

O_3 = Operator(oracleCircuit_3).data
pm(O_3)

# all possible states
state_00 = Kron(state_0, state_0)
state_01 = Kron(state_0, state_1)
state_10 = Kron(state_1, state_0)
state_11 = Kron(state_1, state_1)

print("O|00>", O_3 @ state_00)
print("O|01>", O_3 @ state_01)
print("O|10>", O_3 @ state_10)
print("O|11>", O_3 @ state_11) # flip phase

# H R H, where R = 2|0><0| - I 

def createR_2q():
    circuit = QuantumCircuit(2, 2)
    
    circuit.z(0)
    circuit.z(1)
    circuit.cz(0, 1)
    
    return circuit

R_2q = createR_2q()
R_2q.draw(output='mpl')

R_2q = Operator(R_2q).data
pm(R_2q)

# flip phase(no work on zero states) => Conditional Phase Shift gate
print("RO|00>", R_2q @ O_3 @ state_00) # zero state => would not flip
print("RO|01>", R_2q @ O_3 @ state_01)
print("RO|10>", R_2q @ O_3 @ state_10)
print("RO|11>", R_2q @ O_3 @ state_11)

def createDiffuser_2q():
    circuit = QuantumCircuit(2, 2)
    
    circuit.h(0)
    circuit.h(1)
    circuit = circuit.compose(createR_2q())
    circuit.h(0)
    circuit.h(1)
    circuit.barrier()
    
    return circuit

diffuserCircuit_2q = createDiffuser_2q()
diffuserCircuit_2q.draw(output='mpl')

diff_2q = Operator(diffuserCircuit_2q).data
pm(diff_2q)

DB = 0.5 * state_00 + 0.5 * state_01 + 0.5 * state_10 + 0.5 * state_11
print("DO|DB>", diff_2q @ O_3 @ DB)

def createGroverIteration(oracle, diffuser):
    return oracle.compose(diffuser)

groverIteration_2q = createGroverIteration(createOracle_3(), createDiffuser_2q())
groverIteration_2q.draw(output='mpl')

groverIteration_2q = createGroverIteration(createOracle_3(), createDiffuser_2q())

grover_2q_1 = initCircuit(2).compose(groverIteration_2q.copy())
grover_2q_1.draw(output='mpl')

grover_2q_1.measure([0, 1], [0, 1])

job = execute(grover_2q_1, simulator, shots = 10000)
results = job.result()
counts = results.get_counts(grover_2q_1)
print(counts)

plot_histogram(counts, figsize=(1, 5), color="#CC3333", title="one iteration - find 3")

grover_2q_2 = initCircuit(2).compose(groverIteration_2q.copy()).compose(groverIteration_2q.copy())
grover_2q_2.draw(output='mpl')

grover_2q_2.measure([0, 1], [0, 1])

job = execute(grover_2q_2, simulator, shots = 10000)
results = job.result()
counts = results.get_counts(grover_2q_2)
print(counts)

plot_histogram(counts, figsize=(7, 5), color="#FF9999", title="two iteration - find 3")

# https://www.quantum-inspire.com/kbase/grover-algorithm/

inputCircuit_3q = initCircuit(3)
inputCircuit_3q.draw(output='mpl')

def createOracle_6():
    circuit = QuantumCircuit(3, 3)

    # Oracle for find 6

    # U_f 
    circuit.x(0)
    circuit.h(2)
    circuit.ccx(0, 1, 2)
    circuit.h(2)
    circuit.x(0)
    
    circuit.barrier()
    
    return circuit

oracleCircuit_6 = createOracle_6()
oracleCircuit_6.draw(output='mpl')

# where the entry that correspond to the marked item will have a negative phase
pm(Operator(oracleCircuit_6).data)

Operator(oracleCircuit_6).data

# H R H, where R = 2|0><0| - I 

def createR_3q():
    circuit = QuantumCircuit(3, 3)
    
    circuit.x(2)
    circuit.x(0)
    circuit.x(1)
    circuit.h(2)
    circuit.ccx(0, 1, 2)
    circuit.barrier(0)
    circuit.barrier(1)
    circuit.h(2)
    circuit.x(2)
    circuit.x(0)
    circuit.x(1)
    
    return circuit

R_3q = createR_3q()
R_3q.draw(output='mpl')

pm(Operator(R_3q).data)

print('|000>', Kron(state_0, state_0, state_0))
print('R|000>',Operator(R_3q).data @ Kron(state_0, state_0, state_0) )

print('|010>', Kron(state_0, state_1, state_0))
print('R|010>',Operator(R_3q).data @ Kron(state_0, state_1, state_0) )

print('|110>', Kron(state_1, state_1, state_0))
print('R|110>',Operator(R_3q).data @ Kron(state_1, state_1, state_0) )

def createDiffuser_3q():
    circuit = QuantumCircuit(3, 3)
    
    circuit.h(0)
    circuit.h(1)
    circuit.h(2)
    circuit = circuit.compose(createR_3q())
    circuit.h(0)
    circuit.h(1)
    circuit.h(2)
    circuit.barrier()
    
    return circuit

diffuserCircuit_3q = createDiffuser_3q()
diffuserCircuit_3q.draw(output='mpl')

pm(Operator(diffuserCircuit_3q).data)

groverIteration_3q = createGroverIteration(createOracle_6(), createDiffuser_3q())

groverIteration_3q.draw(output='mpl')

inputCircuit_3q = initCircuit(3)
inputCircuit_3q.draw(output='mpl')

inputCircuit_3q.measure([0, 1, 2], [0, 1, 2])

job = execute(inputCircuit_3q, simulator, shots = 10000)
results = job.result()
counts = results.get_counts(inputCircuit_3q)
print(counts)

plot_histogram(counts, figsize=(10, 5), color="#FFFFCC", title="superposition state - 3 qubits")

grover_3q_1 = initCircuit(3).compose(groverIteration_3q.copy())
grover_3q_1.draw(output='mpl')

grover_3q_1.measure([0, 1, 2], [0, 1, 2])

job = execute(grover_3q_1, simulator, shots = 10000)
results = job.result()
counts = results.get_counts(grover_3q_1)
print(counts)

plot_histogram(counts, figsize=(10, 5), color="#FFBB00", title="one iteration - find 6")

grover_3q_2 = initCircuit(3).compose(groverIteration_3q.copy()).compose(groverIteration_3q.copy())
# grover_3q_2.draw(output='mpl')

grover_3q_2.measure([0, 1, 2], [0, 1, 2])

job = execute(grover_3q_2, simulator, shots = 10000)
results = job.result()
counts = results.get_counts(grover_3q_2)
print(counts)

plot_histogram(counts, figsize=(10, 5), color="#FF8822", title="two iteration - find 6")

grover_3q_3 = initCircuit(3).compose(groverIteration_3q.copy()).compose(groverIteration_3q.copy()).compose(groverIteration_3q.copy())
# grover_3q_3.draw(output='mpl')

grover_3q_3.measure([0, 1, 2], [0, 1, 2])

job = execute(grover_3q_3, simulator, shots = 10000)
results = job.result()
counts = results.get_counts(grover_3q_3)
print(counts)

plot_histogram(counts, figsize=(10, 5), color="#FFFF3F", title="three iteration - find 6")

import math

theta = math.degrees( math.asin(1 / math.sqrt(8)) )
theta

k_tilde = 180 / (4 * theta) - 0.5
k_tilde

for k in range(1, 4):
    print(f'k = {k}: { math.sin(math.radians( (2 * k+1) * theta) ) ** 2 }')

# H R H, where R = 2|0><0| - I 
# A R A^(-1)

def initAACircuit(n):
    circuit = QuantumCircuit(n, n)
    
    circuit.x(0)
    circuit.h(1)
    circuit.z(1)
    circuit.h(0)
    circuit.barrier()
    
    return circuit

inputAACircuit = initAACircuit(2)
inputAACircuit.draw(output='mpl')

A = Operator(initAACircuit(2))

np.allclose(A.data @ Dag(A.data), np.eye(4))

inputAACircuit.measure([0, 1], [0, 1])

job = execute(inputAACircuit, simulator, shots = 10000)
results = job.result()
counts = results.get_counts(inputAACircuit)
print(counts)

plot_histogram(counts, figsize=(7, 5), color="#CCCC66", title="measurements result of A")

oracleCircuit_3 = createOracle_3()
oracleCircuit_3.draw(output='mpl')

def createAADiffuser():
    circuit = QuantumCircuit(2, 2)
    
    circuit.append(A, [0, 1])
    circuit = circuit.compose(createR_2q())
    circuit.append(A.conjugate(), [0, 1])
    circuit.barrier()
    
    return circuit

diffuserAACircuit = createAADiffuser()
diffuserAACircuit.draw(output='mpl')

aaGroverIteration = createGroverIteration(createOracle_3(), createAADiffuser())
aaGroverIteration.draw(output='mpl')

incorrectGroverIteration = createGroverIteration(createOracle_3(), createDiffuser_2q())

grover_incorrect_diffuser = initAACircuit(2).compose(incorrectGroverIteration)
grover_incorrect_diffuser.draw(output='mpl')

grover_incorrect_diffuser.measure([0, 1], [0, 1])

job = execute(grover_incorrect_diffuser, simulator, shots = 10000)
results = job.result()
counts = results.get_counts(grover_incorrect_diffuser)
print(counts)

plot_histogram(counts, figsize=(1, 5), color="#99CC33", title="using incorrect diffuser - find 3")

grover_correct_diffuser = initAACircuit(2).compose(aaGroverIteration)
grover_correct_diffuser.draw(output='mpl')

grover_correct_diffuser.measure([0, 1], [0, 1])

job = execute(grover_correct_diffuser, simulator, shots = 10000)
results = job.result()
counts = results.get_counts(grover_correct_diffuser)
print(counts)

plot_histogram(counts, figsize=(1, 5), color="#339933", title="using correct diffuser - find 3")

inputCircuit_2q = initCircuit(2)
inputCircuit_2q.draw(output='mpl')

def createEvenOracle():
    circuit = QuantumCircuit(2, 2)
    
    circuit.x(0)
    circuit.x(1)
    circuit.cz(0, 1)
    circuit.x(1)
    circuit.cz(0, 1)
    circuit.x(0)
    circuit.barrier()
    
    return circuit

evenOracle = createEvenOracle()
evenOracle.draw(output='mpl')

Operator(evenOracle).data # find 0 and 2

diffuserCircuit_2q = createDiffuser_2q()
diffuserCircuit_2q.draw(output='mpl')

grover_even = initCircuit(2).compose(createEvenOracle()).compose(createDiffuser_2q())
grover_even.draw(output='mpl')

grover_even.measure([0, 1], [0, 1])

job = execute(grover_even, simulator, shots = 10000)
results = job.result()
counts = results.get_counts(grover_even)
print(counts)

plot_histogram(counts, figsize=(5, 5), color="#66CCCC", title="origin grover - find evens")

initState = Kron(
    state_0, state_0, state_0, state_0, state_0, 
    state_0, state_0, state_0, state_0, state_0
)

H = np.array([[1,1], [1,-1]]) / (np.sqrt(2))
H_9 = Kron(H, H, H, H, H, H, H, H, H)
H_10 = Kron(H, H, H, H, H, H, H, H, H, H)

def createTargetState(n):
    states = []
    s = bin(n)[2:].rjust(10, "0")
    for i in list(s):
        if i == '0':
            states.append(state_0)
        else:
            states.append(state_1)
    return Kron(*states)

def createTargetStates(p):
    n = math.floor(p * 1024)
    
    state = createTargetState(0)
    for i in range(1, n):
        state = state + createTargetState(i)
        
    return state / (np.sqrt(n))

target_states = [
    [createTargetStates(0.05), 0.05], # 0.05
    [createTargetStates(0.1), 0.1], # 0.1
    [createTargetStates(0.15), 0.15], # 0.15
    [createTargetStates(0.2), 0.2], # 0.2
    [createTargetStates(0.25), 0.25], # 0.25
    [createTargetStates(0.3), 0.3], # 0.3
    [createTargetStates(0.35), 0.35], # 0.35
    [createTargetStates(0.4), 0.4], # 0.4
    [createTargetStates(0.45), 0.45], # 0.45
    [createTargetStates(0.5), 0.5], # 0.5
    [createTargetStates(0.55), 0.55], # 0.55
    [createTargetStates(0.6), 0.6], # 0.6
    [createTargetStates(0.65), 0.65], # 0.65
    [createTargetStates(0.7), 0.7], # 0.7
    [createTargetStates(0.75), 0.75], # 0.75
    [createTargetStates(0.8), 0.8], # 0.8
    [createTargetStates(0.85), 0.85], # 0.85
    [createTargetStates(0.9), 0.9], # 0.9
    [createTargetStates(0.95), 0.95], # 0.95
    [createTargetStates(1), 1] # 1
]

def createOriginalOracle(target_state):
    N = 1024
    O = np.eye(N) - 2 * ( Dag(target_state).reshape(N, 1) @ target_state.reshape(1, N) )
    
    return O

def createOriginalDiffuser():
    N = 1024

    R = 2 * ( 
        Dag(initState).reshape(N, 1) @ 
        initState.reshape(1, N) 
    ) - np.eye(N)
    
    return  H_10 @ R @ Dag(H_10)

def createHalfPIOracle(target_state):
    N = 1024
    alpha = np.pi / 2
    alpha_phase = np.exp(alpha * 1j)

    # Oracle = I - (1 - e^(αi))|target><target|

    O = np.eye(N) - (1 - alpha_phase) * \
        ( Dag(target_state).reshape(N, 1) @ target_state.reshape(1, N) )
    
    return O

def createHalfPIDiffuser():
    N = 1024
    beta = -np.pi / 2
    beta_phase = np.exp(beta * 1j)

    # R = (1 - e^(βi))|0><0| + e^(βi)I
    R = (1 - beta_phase) * ( 
        Dag(initState).reshape(N, 1) @ 
        initState.reshape(1, N) 
    ) + beta_phase * np.eye(N)
    
    return  H_10 @ R @ Dag(H_10)

def createOneTenthPIOracle(target_state):
    N = 1024
    alpha = np.pi / 10
    alpha_phase = np.exp(alpha * 1j)

    # Oracle = I - (1 - e^(αi))|target><target|

    O = np.eye(N) - (1 - alpha_phase) * \
        ( Dag(target_state).reshape(N, 1) @ target_state.reshape(1, N) )
    
    return O

def createOneTenthPIDiffuser():
    N = 1024
    beta = -np.pi / 10
    beta_phase = np.exp(beta * 1j)

    # R = (1 - e^(βi))|0><0| + e^(βi)I
    R = (1 - beta_phase) * ( 
        Dag(initState).reshape(N, 1) @ 
        initState.reshape(1, N) 
    ) + beta_phase * np.eye(N)
    
    return  H_10 @ R @ Dag(H_10)

def createYouneOracle(n):
    O = np.eye(1024)
    for i in range(n):
        O[i][i] = 0
        O[i][i + 512] = 1
        O[i + 512][i + 512] = 0
        O[i + 512][i] = 1
    
    return O

def createYouneDiffuser():
    N = 512

    R = (2 * Dag(initState).reshape(N * 2, 1) @ 
        initState.reshape(1, N * 2) - np.eye(N * 2))
    
    return  Kron(np.eye(2), H_9) @ R @ Dag(Kron(np.eye(2), H_9))

measurements = getMeasurements(10)

# test: one grover iteration
x_o = []
y_o = []
for target_state, p in target_states:
    n = math.floor(p * 1024)

    O = createOriginalOracle(target_state)
    diff = createOriginalDiffuser()
    final_state = diff @ O @ H_10 @ initState

    probability = 0.0
    for i in range(1024):
        if i < n:
            probability = probability + measure(measurements[i], final_state)
    
#     print(f'p={p}, probability={probability}')
    x_o.append(p)
    y_o.append(probability)

plt.title("Original Grover - apply one iterator")
plt.scatter(x_o, y_o, color="#66CCCC")

# test: best probability
x_o_best = []
y_o_best = []
for target_state, p in target_states:
    n = math.floor(p * 1024)

    O = createOriginalOracle(target_state)
    diff = createOriginalDiffuser()
    final_state = H_10 @ initState

    theta = math.degrees( math.asin(math.sqrt(p)) )
    k_tilde = math.floor( 180 / (4 * theta) - 0.5 )
    for i in range(k_tilde):
        final_state = diff @ O @ final_state

    probability = 0.0
    for i in range(1024):
        if i < n:
            probability = probability + measure(measurements[i], final_state)
    
    final_state = diff @ O @ final_state
    probability2 = 0.0
    for i in range(1024):
        if i < n:
            probability2 = probability2 + measure(measurements[i], final_state)
            
#     print(f'p={p}, probability={probability}')
    x_o_best.append(p)
    y_o_best.append(max(probability, probability2))

plt.title("Original Grover - best Probability")
x = np.arange(-0,1,0.01)
y = 1 - x
plt.plot(x, y, ls='--', color="#FFCC00")

plt.ylim(0, 1.05)
plt.scatter(x_o_best, y_o_best, color="#66CCCC")

# test: one grover iteration
x_half_pi = []
y_half_pi= []
for target_state, p in target_states:
    n = math.floor(p * 1024)

    O = createHalfPIOracle(target_state)
    diff = createHalfPIDiffuser()
    final_state = diff @ O @ H_10 @ initState

    probability = 0.0
    for i in range(1024):
        if i < n:
            probability = probability + measure(measurements[i], final_state)
    
#     print(f'p={p}, probability={probability}')
    x_half_pi.append(p)
    y_half_pi.append(probability)

plt.title("pi/2 Grover - apply one iterator")
plt.scatter(x_half_pi, y_half_pi, color="#66CCCC")

# one grover iteration - pi/2 vs origin

fig, ax = plt.subplots()
plt.title("Apply one iteration")
plt.axvline(1/3, linestyle='--', color="#FFCCCC", label="t/N = 1/3")
plt.axhline(0.925, linestyle='--', color="#FFCC00", label="P=0.925")
plt.scatter(x_o, y_o, color="#66CCCC", label="Original Grover")
plt.scatter(x_half_pi, y_half_pi, color="#FF6666", label="pi/2 phase based Grover")
ax.legend()

# test: one grover iteration
x_ot_pi = []
y_ot_pi= []
for target_state, p in target_states:
    n = math.floor(p * 1024)

    O = createOneTenthPIOracle(target_state)
    diff = createOneTenthPIDiffuser()
    final_state = diff @ O @ H_10 @ initState

    probability = 0.0
    for i in range(1024):
        if i < n:
            probability = probability + measure(measurements[i], final_state)
    
#     print(f'p={p}, probability={probability}')
    x_ot_pi.append(p)
    y_ot_pi.append(probability)

plt.title("0.1pi Grover - apply one iteration")
plt.scatter(x_ot_pi, y_ot_pi, color="#66CCCC")

# test: best probability
x_ot_best = []
y_ot_best = []
for target_state, p in target_states:
    n = math.floor(p * 1024)

    O = createOneTenthPIOracle(target_state)
    diff = createOneTenthPIDiffuser()
    final_state = H_10 @ initState

    p_best = 0.0
    k_tilde = math.floor( 5 * np.sqrt(1024) )
    for i in range(k_tilde):
        final_state = diff @ O @ final_state

        probability = 0.0
        for i in range(1024):
            if i < n:
                probability = probability + measure(measurements[i], final_state)
                
        p_best = max(p_best, probability)
    
#     print(f'p={p}, probability={probability}')
    x_ot_best.append(p)
    y_ot_best.append(p_best)

plt.title("0.1pi Grover - best Probability")
plt.ylim(0, 1)
plt.scatter(x_ot_best, y_ot_best, color="#66CCCC")

# best probability - 0.1pi vs origin

fig, ax = plt.subplots()
plt.title("Best Probability")
x = np.arange(-0,1,0.01)
y = 1 - x
plt.plot(x, y, ls='--', color="#FFCC00")

plt.ylim(0, 1.05)
plt.scatter(x_o_best, y_o_best, color="#66CCCC", label="Original Grover")
plt.scatter(x_ot_best, y_ot_best, color="#FF6666", label="0.1pi phase based Grover")
ax.legend()

# test: best probability
x_y_best = []
y_y_best = []
for _, p in target_states:
    n = math.floor(p * 512)

    O = createYouneOracle(n)
    diff = createYouneDiffuser()
    final_state = Kron( np.eye(2), H_9) @ initState
    
    p_best = 0.0
    k_tilde = math.floor(np.pi / (2 * np.sqrt(2)) * np.sqrt(512 / n))
    for i in range(k_tilde):
        final_state = diff @ O @ final_state
        
    probability = 0.0
    for i in range(512):
        if i < n:
            probability = probability + measure(measurements[i], final_state) + measure(measurements[i + 512], final_state)

    #     print(f'p={p}, probability={probability}')
    x_y_best.append(p)
    y_y_best.append(probability)

fig, ax = plt.subplots()
plt.title("Best Probability")

x = np.arange(-0, 1, 0.01)
y = 1 - x
plt.plot(x, y, ls='--', color="#FFCC00")

plt.axhline(0.8472, linestyle='--', color="#99CCFF", label="P=0.8472")
plt.axvline(0.308, linestyle='--', color="#FFCCCC", label="P=0.308")

plt.ylim(0, 1.05)
plt.scatter(x_y_best, y_y_best, color="#FF6666", label="Youne")
plt.scatter(x_o_best, y_o_best, color="#66CCCC", label="Original Grover")
ax.legend()

def createYouneInitCircuit():
    circuit = QuantumCircuit(3, 2)
    
    circuit.h(0)
    circuit.h(1)
    circuit.barrier()
    
    return circuit

youneInitCircuit = createYouneInitCircuit()
youneInitCircuit.draw(output='mpl')

def createYouneOracle():
    circuit = QuantumCircuit(3, 2)
    
    circuit.x(0)
    circuit.x(1)
    circuit.ccx(0, 1, 2)
    circuit.x(0)
    circuit.x(1)
    circuit.barrier()
    
    circuit.x(0)
    circuit.ccx(0, 1, 2)
    circuit.x(0)
    circuit.barrier()
    
    return circuit

youneOracle = createYouneOracle()
youneOracle.draw(output='mpl')

def createYouneDiffuser():
    circuit = QuantumCircuit(3, 2)
    
    circuit.h(0)
    circuit.h(1)
    circuit.barrier(2)
    
    circuit.x(2)
    circuit.x(0)
    circuit.x(1)
    circuit.h(2)
    circuit.ccx(0, 1, 2)
    circuit.barrier(0)
    circuit.barrier(1)
    circuit.h(2)
    circuit.x(2)
    circuit.x(0)
    circuit.x(1)
    
    circuit.h(0)
    circuit.h(1)
#     circuit.h(2)
    circuit.barrier()
    
    return circuit

youneDiffuser = createYouneDiffuser()
youneDiffuser.draw(output='mpl')

youneGroverIteration = createYouneOracle().compose(createYouneDiffuser())

grover_youne = createYouneInitCircuit().compose(youneGroverIteration.copy())
grover_youne.draw(output='mpl')

grover_youne.measure([0, 1], [0, 1])

job = execute(grover_youne, simulator, shots = 10000)
results = job.result()
counts = results.get_counts(grover_youne)
print(counts)

plot_histogram(counts, figsize=(4, 5), color="#0099CC", title="youne grover - find evens")

def createYouneOracle2():
    circuit = QuantumCircuit(3, 2)
    
    circuit.x(0)
    circuit.cx(0, 2)
    circuit.x(0)
    circuit.barrier()
    
    return circuit

youneOracle2 = createYouneOracle2()
youneOracle2.draw(output='mpl')

# 0 => 4 2 => 6
Operator(youneOracle2).data

youneGroverIteration2 = createYouneOracle2().compose(createYouneDiffuser())

grover_youne2 = createYouneInitCircuit().compose(youneGroverIteration2.copy())
grover_youne2.draw(output='mpl')

grover_youne2.measure([0, 1], [0, 1])

job = execute(grover_youne2, simulator, shots = 10000)
results = job.result()
counts = results.get_counts(grover_youne2)
print(counts)

plot_histogram(counts, figsize=(4, 5), color="#66CCFF", title="youne grover(using another Oracle) - find evens")

# https://towardsdatascience.com/behind-oracles-grovers-algorithm-amplitude-amplification-46b928b46f1e

def createSATInitCircuit(var, clause):
    # input + workspace(clause) + checker(1)
    circuit = QuantumCircuit(var + clause + 1, var)
    
    for i in range(var):
        circuit.h(i)
    circuit.barrier()
    
    return circuit

SATInitCircuit = createSATInitCircuit(4, 3)
SATInitCircuit.draw(output='mpl')

def createSATOracle(var, clause):
    circuit = QuantumCircuit(var + clause + 1, var)
    
    # (a ∧ b ∧ ¬c)
    circuit.x(2)
    circuit.mcx([0, 1, 2], 4)
    circuit.x(2)
    circuit.barrier()
    
    # (¬b ∧d )
    circuit.x(1)
    circuit.mcx([1, 3], 5)
    circuit.x(1)
    circuit.barrier()
    
    # ¬d
    circuit.x(3)
    circuit.cx(3, 6)
    circuit.x(3)
    circuit.barrier()
    
    # (a ∧ b ∧ ¬c) ∧ ¬(¬b ∧d ) ∧ ¬d
    circuit.x(5)
    circuit.mcx([4, 5, 6], 7)
    circuit.x(5)
    circuit.barrier()
    
    # uncomputation
    
    # ¬d
    circuit.x(3)
    circuit.cx(3, 6)
    circuit.x(3)
    circuit.barrier()
    
    # (¬b ∧d )
    circuit.x(1)
    circuit.mcx([1, 3], 5)
    circuit.x(1)
    circuit.barrier()
    
    # (a ∧ b ∧ ¬c)
    circuit.x(2)
    circuit.mcx([0, 1, 2], 4)
    circuit.x(2)
    circuit.barrier()
    
    return circuit

SATOracleCircuit = createSATOracle(4, 3)
SATOracleCircuit.draw(output='mpl')

def createSATDiffuser(var, clause):
    circuit = QuantumCircuit(var + clause + 1, var)
    
    circuit.h(0)
    circuit.h(1)
    circuit.h(2)
    circuit.h(3)
    circuit.x(0)
    circuit.x(1)
    circuit.x(2)
    circuit.x(3)
    circuit.mcx([0, 1, 2, 3], var + clause)
    circuit.x(0)
    circuit.x(1)
    circuit.x(2)
    circuit.x(3)
    circuit.h(0)
    circuit.h(1)
    circuit.h(2)
    circuit.h(3)
    
    return circuit

SATDiffuserCircuit = createSATDiffuser(4, 3)
SATDiffuserCircuit.draw(output='mpl')

SATCircuit = createSATInitCircuit(4, 3).compose(createSATOracle(4, 3)).compose(createSATDiffuser(4, 3))
SATCircuit.draw(output='mpl')

SATCircuit.measure([0, 1, 2, 3], [0, 1, 2, 3])

job = execute(SATCircuit, simulator, shots = 10000)
results = job.result()
counts = results.get_counts(SATCircuit)
print(counts)

plot_histogram(counts, figsize=(12, 5), color="#CCCCFF", title="SAT - solving (a ∧ b ∧ ¬c) ∧ ¬(¬b ∧d ) ∧ ¬d, k = 1")

SATIteration = Operator(createSATOracle(4, 3).compose(createSATDiffuser(4, 3)))

SATCircuit2 = createSATInitCircuit(4, 3)
SATCircuit2.append(SATIteration, list(range(8)))
SATCircuit2.append(SATIteration, list(range(8)))
SATCircuit2.append(SATIteration, list(range(8)))
SATCircuit2.draw(output='mpl')

SATCircuit2.measure([0, 1, 2, 3], [0, 1, 2, 3])

job = execute(SATCircuit2, simulator, shots = 10000)
results = job.result()
counts = results.get_counts(SATCircuit2)
print(counts)

plot_histogram(counts, figsize=(12, 5), color="#6666FF", title="SAT - solving (a ∧ b ∧ ¬c) ∧ ¬(¬b ∧d ) ∧ ¬d, k = 3")

import math

math.degrees( math.asin(1 / 4))

for k in range(1, 5):
    print(f'k = {k}: { math.sin(math.radians( (2 * k+1) * 14.47) ) ** 2 }')

# https://arxiv.org/pdf/quant-ph/9605034.pdf

inputCircuit_8q = initCircuit(8)
inputCircuit_8q.draw(output='mpl')

def createOracle_255():
    circuit = QuantumCircuit(8, 8)
    
    circuit.h(7)
    circuit.mcx(list(range(7)), 7)
    circuit.h(7)
    circuit.barrier()
    
    return circuit

oracleCircuit_255 = createOracle_255()
oracleCircuit_255.draw(output='mpl')

Operator(oracleCircuit_255).data

for index, row in enumerate(Operator(oracleCircuit_255).data):
    if abs(row[index] - -1) < 1e-6:
        print(index)

def createR_8q():
    circuit = QuantumCircuit(8, 8)
    
    circuit.x(7)
    circuit.x(6)
    circuit.x(5)
    circuit.x(4)
    circuit.x(3)
    circuit.x(2)
    circuit.x(0)
    circuit.x(1)
    circuit.h(7)
    circuit.mcx(list(range(7)), 7)
    circuit.barrier(0)
    circuit.barrier(1)
    circuit.barrier(2)
    circuit.barrier(3)
    circuit.barrier(4)
    circuit.barrier(5)
    circuit.barrier(6)
    circuit.h(7)
    circuit.x(2)
    circuit.x(0)
    circuit.x(1)
    circuit.x(5)
    circuit.x(4)
    circuit.x(3)
    circuit.x(7)
    circuit.x(6)
    
    return circuit

R_8q = createR_8q()
R_8q.draw(output='mpl')

def createDiffuser_8q():
    circuit = QuantumCircuit(8, 8)
    
    circuit.h(0)
    circuit.h(1)
    circuit.h(2)
    circuit.h(3)
    circuit.h(4)
    circuit.h(5)
    circuit.h(6)
    circuit.h(7)
    circuit = circuit.compose(createR_8q())
    circuit.h(0)
    circuit.h(1)
    circuit.h(2)
    circuit.h(3)
    circuit.h(4)
    circuit.h(5)
    circuit.h(6)
    circuit.h(7)
    circuit.barrier()
    
    return circuit

diffuserCircuit_8q = createDiffuser_8q()
diffuserCircuit_8q.draw(output='mpl')

groverIteration_8q = createGroverIteration(createOracle_255(), createDiffuser_8q())
groverIteration_8q.draw(output='mpl')

fullCircuit_8q = initCircuit(8).compose(groverIteration_8q)
fullCircuit_8q.draw(output='mpl')

degree = math.degrees( math.asin(1 / np.sqrt(256)) )
degree

for k in range(1, 15):
    print(f'k = {k}: { math.sin(math.radians( (2 * k+1) * degree) ) ** 2 }')

def groverMeasurements(k):
    circuit = initCircuit(8)
    
    for i in range(k):
        circuit = circuit.compose(groverIteration_8q.copy())
    
    circuit.measure(list(range(8)), list(range(8)))
    
    job = execute(circuit, simulator, shots = 1)
    results = job.result()
    counts = results.get_counts(circuit)
    
    return list(counts)[0]

from random import randrange

def solutionFinding(limit):
    m = 1
    lam = 6 / 5
    count = 0 # iteration times
    
    while count <= limit * (9/4):
        if m > 1:
            j = randrange(0, math.floor(m))
            if j != 0:
                count += j
                answer = groverMeasurements(j)

                if answer == '11111111':
#                     print(f'Find! - total times: {count}')
                    return count, True
        m = min(lam * m, limit)
        
#     print(f'Not find! - total times: {count}')
    return count, False

times = 0
best = 1000
worst = 0
success = 0
for i in range(1000):
    time, flag = solutionFinding(np.sqrt(256)) # √N
    times += time
    best = min(best, time)
    worst = max(worst, time)
    if flag:
        success += 1
print(f'total test times: {1000}')
print('--------------------------')
print(f'Average run times: {times / 1000}')
print(f'best run times: {best}')
print(f'worst run times: {worst}')
print(f'Success rate: {success / 1000 * 100:.2f}%')

