# imports
import math

# importing Qiskit
from qiskit import QuantumCircuit, Aer, assemble, transpile
from qiskit.visualization import plot_histogram

f = open("cnf_test1.txt")
info = f.readline().split(' ')
literals = int(info[0])
clausulas = int(info[1])

lines = [line.rstrip('\n') for line in f]

formula = []
for i in lines:
    l = i.split(' ')
    formula.append([int(l[0]),int(l[1]),int(l[2])])

n = literals + clausulas + 2 # Number of qubits in the circuit

def initialize_grover(qc, qubits):
    # Superposition
    for i in range(literals+1):
        qc.h(qubits[i])
        
    # State |-⟩
    qc.x(n-1)
    qc.h(n-1)
    
    return qc

qc = QuantumCircuit(n)

# Creation of each clause.
contador = 0
for c in formula:
    for i in range(3):
        if c[i] > 0:
            qc.x(abs(c[i])-1)
            
    qc.mct(list([abs(c[0])-1,abs(c[1])-1,abs(c[2])-1]), literals+1+contador)  # multi-controlled-toffoli
    
    for i in range(3):
        if c[2-i] > 0:
            qc.x(abs(c[2-i])-1)
    
    contador+=1


# Checking of each clause
qc.x(literals) # Augmented space
for q in range(literals+1,n-1):
    qc.x(q)
    
qc.mct(list(list([literals, *range(literals+1,n-1)])), n-1)

for q in range(literals+1,n-1):
    qc.x(q)
qc.x(literals) # Augmented space


# Creation of each clause (symmetry).    
contador -= 1
for c in reversed(formula):
    for i in range(3):
        if c[i] > 0:
            qc.x(abs(c[i])-1)
            
    qc.mct(list([abs(c[0])-1,abs(c[1])-1,abs(c[2])-1]), literals+1+contador)
    
    for i in range(3):
        if c[2-i] > 0:
            qc.x(abs(c[2-i])-1)
    
    contador-=1

oracle = qc.to_gate()
oracle.name = "Oracle"
qc.draw()

qc = QuantumCircuit(n)

for qubit in range(literals+1):
        qc.h(qubit)
        
for qubit in range(literals+1):
        qc.x(qubit)
        
# Equivalent Diffuser. It gives us the opposite result in Quantum counting, therefore M will be the number of no solutions at the end.
qc.mct(list([i for i in range(literals+1)]), n-1)

########## Original Diffuser ##############################################################
#qc.mct(list([i for i in range(literals+1)]), literals+1)
#qc.x(literals+1)
#qc.cx(literals+1,n-1)
#qc.x(literals+1)
#qc.mct(list([i for i in range(literals+1)]), literals+1)
###########################################################################################

for qubit in range(literals+1):
        qc.x(qubit)
        
for qubit in range(literals+1):
        qc.h(qubit)

inversionMean = qc.to_gate()
inversionMean.name = "inversionMean"
qc.draw()

grover_iteration = QuantumCircuit(n)

grover_iteration.append(oracle, list(range(n)))
grover_iteration.append(inversionMean, list(range(n)))

gi = grover_iteration.to_gate()
gi.label = "GroverI"
cgi = gi.control()

def qft(n):
    circuit = QuantumCircuit(n)
    
    def swap_registers(circuit, n):
        for qubit in range(n//2):
            circuit.swap(qubit, n-qubit-1)
        return circuit
    
    def qft_rotations(circuit, n):
        #Performs qft on the first n qubits in circuit (without swaps)
        if n == 0:
            return circuit
        n -= 1
        circuit.h(n)
        for qubit in range(n):
            circuit.cp(math.pi/2**(n-qubit), qubit, n)
        qft_rotations(circuit, n)
    
    qft_rotations(circuit, n)
    swap_registers(circuit, n)  # Comment this line if we use the 2 possible order of the QPE registers
    return circuit

t = literals+1  # Number of counting qubits.
qft_dagger = qft(t).to_gate().inverse()
qft_dagger.label = "QFT†"

qpe = QuantumCircuit(n+t, t) # Circuit with n+t qubits and t classical bits

# Superposition
for qubit in range(t+literals+1):
    qpe.h(qubit)

# State |-⟩
qpe.x(n+t-1)
qpe.h(n+t-1)
    
# Grover iterations
iterations = 1
for qubit in range(t):
    for i in range(iterations):
        qpe.append(cgi, [qubit] + [*range(t, n+t)])   # Using this line, uncomment the QFT swap.
        #qpe.append(cgrit, [(t-1)-qubit] + [*range(t, n+t)]) # Using this line, comment the QFT swap.

    iterations *= 2
    
# Application of QFT inverse to all counting qubits
qpe.append(qft_dagger, range(t))

# Measure the counting qubits
qpe.measure(range(t), range(t))

qpe.draw()

# Quantum counting sinulation
aer_sim = Aer.get_backend('aer_simulator')
transpiled_qpe = transpile(qpe, aer_sim)
qobj = assemble(transpiled_qpe)
job = aer_sim.run(qobj)
hist = job.result().get_counts()
plot_histogram(hist)

measured_str = max(hist, key=hist.get)

measured_int = int(measured_str,2)
print(measured_str)
print("Register Output = %i" % measured_int)

theta = (measured_int/(2**t))*math.pi*2
print("Theta = %.5f" % theta)

N = 2**(literals+1)
M = N * (math.sin(theta/2)**2)
M = N-M # Diffuser faster only
print(N,M)
print("No. of Solutions = %.1f" % M)

# Upper bound error for the number of solutions
m = t - 1
err = (math.sqrt(2*M*N) + N/(2**(m+1)))*(2**(-m))
print("Error < %.2f" % err)

rounds = math.pi/4.0*math.sqrt(N/M)
print("Result:",rounds)
print("Iterations R ≤",math.ceil(rounds))

if math.floor(rounds) == 0:
    rounds = round(rounds)
else:
    rounds = math.floor(rounds)

print("Iterations chosen:", rounds)

grover_circuit = QuantumCircuit(n,literals+1)

grover_circuit = initialize_grover(grover_circuit, list(range(n)))
for i in range(rounds):
    grover_circuit.append(oracle, list(range(n)))
    grover_circuit.append(inversionMean, list(range(n)))

for i in range(literals+1):
    grover_circuit.measure([i],[i])

grover_circuit.draw()

# Grover's algorithm simulation
aer_sim = Aer.get_backend('aer_simulator')
transpiled_grover_circuit = transpile(grover_circuit, aer_sim)
qobj = assemble(transpiled_grover_circuit)
results = aer_sim.run(qobj, shots=1024).result()
counts = results.get_counts()
print("\nTotal count:",counts)
plot_histogram(counts,figsize=(12,5))