import qiskit
from qiskit import QuantumCircuit
from qiskit.visualization import plot_histogram
from qiskit.quantum_info import Operator

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')

# 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 }')