#initialization
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

# importing Qiskit
from qiskit import Aer, IBMQ
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister
from qiskit import available_backends, execute, register, get_backend, compile
from qiskit.wrapper.jupyter import *

# import basic plot tools
from qiskit.tools.visualization import plot_histogram

# Load the saved IBMQ accounts
IBMQ.load_accounts()

"""
Explanation: <img src="../../images/qiskit-heading.gif" alt="Note: In order for images to show up in this jupyter notebook you need to select File => Trusted Notebook" width="500 px" align="left">
The Simon Algorithm
The latest version of this notebook is available on https://github.com/qiskit/qiskit-tutorial.
The Simon algorithm is an example that shows a quantum algorithm can solve a problem exponentially efficient than any classical algorithms. Like the Grover search, it depends on the existence of a blackbox (or, oracle) function that returns a predefined output over specific input or query. In the query-complexity setting, one cares only about how many queries are required to solve a specific problem, but does not care how the blackbox is realized. However, in this tutorial we have to implement it using the unit gates available in QISKit, just like we have done with the Grover search.
We first describe the problem addressed by the Simon algorithm, show the steps of the algorithm and the construction of the blackbox function, and present the experimental results on simulators and real devices. 

Contributors
Rudy Raymond
The Problem <a id='introduction'></a>
The Simon algorithm deals with finding a hidden integer $s \in {0,1}^n$ from an oracle $f_s$ that satisfies $f_s(x) = f_s(y)$ if and only if $y = x \oplus s$ for all $x \in {0,1}^n$.  Here, the $\oplus$ is the bitwise XOR operation. Thus, if $s = 0\ldots 0$, i.e., the all-zero bitstring, then $f_s$ is a 1-to-1 (or, permutation) function. Otherwise, if $s \neq 0\ldots 0$, then $f_s$ is a 2-to-1 function. 
The Simon algorithm can find the hidden integer using only $O(n)$ queries to the blackbox function, while any classical algorithms require $\Omega(\sqrt{2^n})$ queries. 
The Algorithm to Find the Hidden Integer
The Simon algorithm finds the hidden integer by combining quantum algorithm with postprocessing on classical computers as below. 


Prepare two quantum registers each of length $n$ that are initialized to all-zero bitstring: the first one as input and the second one as output of the blackbox function.
$$
|0\rangle |0\rangle
$$


Apply Hadamard gates to the first register to create superposition of all possible inputs. 
$$
H^{\otimes n} |0\rangle |0\rangle = \frac{1}{\sqrt{2^n}} \sum_{x=0}^{2^n-1} |x\rangle |0\rangle
$$


Query the blackbox function to obtain the answer to queries on the second register.
$$
\frac{1}{\sqrt{2^n}} \sum_{x=0}^{2^n-1} U_{f_s}|x\rangle |0\rangle = \frac{1}{\sqrt{2^n}} \sum_{x=0}^{2^n-1} |x\rangle |f_s(x)\rangle
$$


Apply Hadamard gates to the first register. 
$$
\frac{1}{\sqrt{2^n}} \sum_{x=0}^{2^n-1} H^{\otimes n}|x\rangle |f_s(x)\rangle = \frac{1}{2^n} \sum_{y=0}^{2^n-1}\sum_{x=0}^{2^n-1} (-1)^{x \cdot y}|y\rangle |f_s(x)\rangle = \frac{1}{2^n} \sum_{y=0}^{2^n-1} |y\rangle \sum_{x=0}^{2^n-1} ( (-1)^{x \cdot y} + (-1)^{(x\oplus s) \cdot y} )  |f_s(x)\rangle
$$
Notice that at the right-hand side of the above equation, because $(-1)^{(x\oplus s) \cdot y} = (-1)^{x\cdot y +  s \cdot y}$ we can conclude that the probability amplitude of the basis state $|y\rangle |f_s(x)\rangle$ is $(-1)^{x\cdot y} (1 + (-1)^{s \cdot y}  )$, which is zero if and only if $s \cdot y = 1$. Thus, measuring the first register will always give $y$ such that $s \cdot y = 0$. Moreover, we can obtain many different $y$'s by repeating Step 1 to 4. 


Repeat Step 1 to 4 for $m$ times to obtain $y_1, y_2, \ldots, y_m$. 


(Classical post-processing) Let $\mathbf{Y}$ be an $m\times n$ matrix whose $i$-th row is $y_i$ in Step 5, and $\vec{s}$ be the column vector whose $j$-th element is the $j$-th bit of $s$. Solve the following system of linear equations to obtain $s$. 
$$
\mathbf{Y} \vec{s} = 0 
$$


The Circuit <a id="circuit"></a>
We now implement the Simon algorithm with Qiskit by first preparing the environment.
End of explanation
"""


s = "010101" # the hidden bitstring

assert 1 < len(s) < 20, "The length of s must be between 2 and 19"

for c in s:
    assert c == "0" or c == "1", "s must be a bitstring of '0' and '1'"
    
n = len(s) #the length of the bitstring

"""
Explanation: We then set the hidden bitstring $s$ that will be used to construct the circuit of the blackbox function (whose details will be given later). The number of qubits used in the experiment is twice the length of the bitstring $s$.
End of explanation
"""


# Step 1
# Creating registers
# qubits for querying the oracle and recording its output
qr = QuantumRegister(2*n)
# for recording the measurement on the first register of qr
cr = ClassicalRegister(n)

circuitName = "Simon"
simonCircuit = QuantumCircuit(qr, cr)

# Step 2
# Apply Hadamard gates before querying the oracle
for i in range(n):
    simonCircuit.h(qr[i])

# Apply barrier to mark the beginning of the blackbox function
simonCircuit.barrier()

"""
Explanation: We then use Qiskit to create the circuit of the Simon algorithm prior the querying the blackbox function.
End of explanation
"""


# Step 3 query the blackbox function
# copy the content of the first register to the second register
for i in range(n):
    simonCircuit.cx(qr[i], qr[n+i])

# get the least index j such that s_j is "1"
j = -1
for i, c in enumerate(s):
    if c == "1":
        j = i
        break

# Creating 1-to-1 or 2-to-1 mapping with the j-th qubit of x as control to XOR the second register with s
for i, c in enumerate(s):
    if c == "1" and j >= 0:
        simonCircuit.cx(qr[j], qr[n+i]) #the i-th qubit is flipped if s_i is 1

# get random permutation of n qubits
perm = list(np.random.permutation(n))

#initial position
init = list(range(n))
i = 0
while i < n:
    if init[i] != perm[i]:
        k = perm.index(init[i])
        simonCircuit.swap(qr[n+i], qr[n+k]) #swap qubits
        init[i], init[k] = init[k], init[i] #marked swapped qubits
    else:
        i += 1
        
# randomly flip the qubit
for i in range(n):
    if np.random.random() > 0.5:
        simonCircuit.x(qr[n+i])

# Apply the barrier to mark the end of the blackbox function
simonCircuit.barrier()

"""
Explanation: Constructing a Circuit for the Blackbox Function
We now details the construction of 1-to-1 and 2-to-1 permutation circuit of the blackbox function. Let us assume the blackbox function receive $|x\rangle|0\rangle$ as input. With regards to a predetermined $s$, the blackbox function writes its output to the second register so that it transforms the input to $|x\rangle|f_s(x)\rangle$ such that $f(x) = f(x\oplus s)$ for all $x \in {0,1}^n$. 
Such a blackbox function can be realized by the following procedures.


Copy the content of the first register to the second register.
$$
|x\rangle|0\rangle \rightarrow |x\rangle|x\rangle
$$


(Creating 1-to-1 or 2-to-1 mapping) If $s$ is not all-zero, then there is the least index $j$ so that $s_j = 1$. If $x_j = 0$, then XOR the second register with $s$. Otherwise, do not change the second register.
$$
|x\rangle|x\rangle \rightarrow |x\rangle|x \oplus s\rangle~\mbox{if}~x_j = 0~\mbox{for the least index j}
$$


(Creating random permutation) Randomly permute and flip the qubits of the second register.
$$
|x\rangle|y\rangle \rightarrow |x\rangle|f_s(y)\rangle
$$


Below is the circuit of the blacbox function based on the above procedures.
End of explanation
"""


# Step 4 apply Hadamard gates to the first register
for i in range(n):
    simonCircuit.h(qr[i])
    
# Step 5 perform measurement on the first register
for i in range(n):
    simonCircuit.measure(qr[i], cr[i])

#draw the circuit
from qiskit.tools.visualization import circuit_drawer
circuit_drawer(simonCircuit)

"""
Explanation: Now we can continue with the steps of the Simon algorithm: applying the Hadamard gates to the first register and measure.
End of explanation
"""


# use local simulator
backend = Aer.get_backend("qasm_simulator")
# the number of shots is twice the length of the bitstring
shots = 2*n
job = execute(simonCircuit, backend=backend, shots=shots)
answer = job.result().get_counts()

plot_histogram(answer)

"""
Explanation: Experimenting with Simulators
We show the experiments of finding the hidden integer with simulators.
End of explanation
"""


# Post-processing step
# Constructing the system of linear equations Y s = 0
# By k[::-1], we reverse the order of the bitstring
lAnswer = [ (k[::-1],v) for k,v in answer.items() if k != "0"*n  ] #excluding the trivial all-zero
#Sort the basis by their probabilities
lAnswer.sort(key = lambda x: x[1], reverse=True)

Y = []
for k, v in lAnswer:
    Y.append( [ int(c) for c in k ] )


#import tools from sympy
from sympy import Matrix, pprint, MatrixSymbol, expand, mod_inverse


Y = Matrix(Y)

#pprint(Y)

#Perform Gaussian elimination on Y
Y_transformed = Y.rref(iszerofunc=lambda x: x % 2==0) # linear algebra on GF(2) 

#to convert rational and negatives in rref of linear algebra on GF(2)
def mod(x,modulus):
    numer, denom = x.as_numer_denom()
    return numer*mod_inverse(denom,modulus) % modulus

Y_new = Y_transformed[0].applyfunc(lambda x: mod(x,2)) #must takecare of negatives and fractional values
#pprint(Y_new)

print("The hidden bistring s[ 0 ], s[ 1 ]....s[",n-1,"] is the one satisfying the following system of linear equations:")
rows, cols = Y_new.shape
for r in range(rows):
    Yr = [ "s[ "+str(i)+" ]" for i, v in enumerate(list(Y_new[r,:])) if v == 1 ]
    if len(Yr) > 0:
        tStr = " + ".join(Yr)
        print(tStr, "= 0")


"""
Explanation: We can see that the results of the measurements are the basis whose inner product with the hidden string $s$ are zero. 
(Notice that the basis on the label of the x-axis in the above plot are numbered from right to left instead of from left to right that we used for $s$)
Gathering the measurement results, we proceed to post-processing with computations that can be done on classical computers. 
Post Processing with Gaussian Elimination
The post processing is done with Gaussian elimination to solve the system of linear equations to determine $s$.
End of explanation
"""


Aer.backends()

IBMQ.backends()

%%qiskit_job_status
#Use one of the available backends 
backend = IBMQ.get_backend("ibmq_16_melbourne")

# show the status of the backend
print("Status of", backend, "is", get_backend(backend).status())

shots = 10*n #run more experiments to be certain
max_credits = 3          # Maximum number of credits to spend on executions. 

job_exp = execute(simonCircuit, backend=backend, shots=shots, max_credits=max_credits)

results = job_exp.result()
answer = results.get_counts(simonCircuit)

plot_histogram(answer, options={'width': 12, 'height': 6})

# Post-processing step
# Constructing the system of linear equations Y s = 0
# By k[::-1], we reverse the order of the bitstring
lAnswer = [ (k[::-1][:n],v) for k,v in answer.items() ] #excluding the qubits that are not part of the inputs
#Sort the basis by their probabilities
lAnswer.sort(key = lambda x: x[1], reverse=True)

Y = []
for k, v in lAnswer:
    Y.append( [ int(c) for c in k ] )

Y = Matrix(Y)

#Perform Gaussian elimination on Y
Y_transformed = Y.rref(iszerofunc=lambda x: x % 2==0) # linear algebra on GF(2)
Y_new = Y_transformed[0].applyfunc(lambda x: mod(x,2)) #must takecare of negatives and fractional values
#pprint(Y_new)

print("The hidden bistring s[ 0 ], s[ 1 ]....s[",n-1,"] is the one satisfying the following system of linear equations:")
rows, cols = Y_new.shape
for r in range(rows):
    Yr = [ "s[ "+str(i)+" ]" for i, v in enumerate(list(Y_new[r,:])) if v == 1 ]
    if len(Yr) > 0:
        tStr = " + ".join(Yr)
        print(tStr, "= 0")

"""
Explanation: As seen above, the system of linear equations is satisfied by the hidden integer $s$. Notice that there can be more than one solutions to the system. In fact, all-zero bitstring is a trivial solution to the system of linear equations. But by having more samples one can narrow down the candidates of the solution, and then test the solution by querying the blackbock in the classical manner.
Experimenting with Real Devices
We see how one can still find out the hidden integer by running the Simon algorithm on real devices. Due to imperfect quantum computers, obtaining the conclusion is not as easy as done with the simulator of perfect quantum computers.
We first register the token and query the available backends.
End of explanation
"""


