!pip install qiskit && pip install pylatexenc && pip install qiskit-aer && pip install qiskit-ibmq-provider

# Useful imports
import numpy as np
import warnings

# Qiskit components
from qiskit import IBMQ
from qiskit import Aer, ClassicalRegister, QuantumRegister, QuantumCircuit, execute
from qiskit.quantum_info.operators import Operator
from qiskit.quantum_info import random_unitary
from qiskit.tools.monitor import backend_overview, backend_monitor

# Visualization
from qiskit.visualization import plot_histogram, plot_state_city
import matplotlib.pyplot as plt
import matplotlib.patches as patch
from matplotlib.path import Path
from matplotlib.animation import FuncAnimation
from IPython.display import Image

# Filter all warnings from now on
warnings.filterwarnings("ignore")

# Define the states
states = ['|00>', '|01>', '|11>', '|10>']
# Define the coordinates of the square's vertices
x = [0, 1, 1, 0, 0]
y = [0, 0, 1, 1, 0]
# Create a plot and set axis limits
plt.figure(figsize=(5, 5))
plt.plot(x, y, 'k-')
plt.scatter(x, y, c='black', s=100)
# Add text labels with variable offset
label_offset = 0.05
for i, state in enumerate(states):
    if i in [0, 1]:
        plt.text(x[i], y[i] - (1.8 * label_offset), state, fontsize=12, ha='center')
    elif i in [2, 3]:
        plt.text(x[i], y[i] + label_offset, state, fontsize=12, ha='center')
# Set axis limits and labels
plt.xlim(-0.25, 1.25)
plt.ylim(-0.25, 1.25)
# Add a title
plt.title('QRW - 2D hypercube')
# Show the plot
plt.grid(False)
plt.show()

# 1 flip of the Grover coin
q = QuantumRegister(6, name='q') # 6 quantum bits, 2 for the coin, 4 for the hypercube
q_circ = QuantumCircuit(q)
for i in range(1):
    q_circ.h([4, 5])
    q_circ.z([4, 5])
    q_circ.cz(4, 5)
    q_circ.h([4, 5])
q_circ.draw('mpl')

# 1 step in the process - a Grover coin toss and a shift operation
q = QuantumRegister(6, name='q') # same setup as before
q_circ = QuantumCircuit(q)
for i in range(1): # we are simulating one shift operation
    q_circ.h([4, 5])
    q_circ.z([4, 5])
    q_circ.cz(4, 5)
    q_circ.h([4, 5])
for i in range(0, 4): # here starts the shift operator, switching one of the four bits according to the rules mentioned above
    q_circ.x(4)
    if i % 2 == 0:
      q_circ.x(5)
    q_circ.ccx(4, 5, i)
q_circ.draw('mpl')

# In this algorithm we have a dimension of n for the state space, then we have log_2(n) states in the coin space --> define the dimensions of the quantum circuit
class QRW_Reg:
    """
    QRW_Reg initializes a Quantum Random Walk registry --> it inputs the variable dim which is the number of qubits in the state space
    dim: number of qubits
    c: the number of coin qubits
    """
    def __init__(self, dim):
        c = (np.log2(dim))
        if c.is_integer() != 1:
            raise TypeError("The number of qubits n have to satisfy the condition that log_2(n) is an integer!")
        self.c = int(c)
        self.dim = dim
        n = int(dim + c)
        self.n = n
        q_r = QuantumRegister(n, 'q')
        cl_r = ClassicalRegister(n, "c")
        self.qr = q_r
        self.cl = cl_r
        q_circuit = QuantumCircuit(self.qr, self.cl)
        self.circuit = q_circuit

# Define the QRW search algorithm as a subclass of the QRW registry
class QRW_Search(QRW_Reg):
    """
    QRW_Search is an implementation of the quantum random walk search algorithm
    based on the paper by N. Shenvi et alter (see References paragraph for more details).
    This class ineherits the QRW_Reg class
    dim: number of qubits in the state space
    state_search: the state that the algorithm searches for
    """
    def __init__(self, dim, state_search):
        QRW_Reg.__init__(self, dim)
        self.circuit = QRW_Reg(dim).circuit
        quantum_reg = QRW_Reg(dim).qr
        classical_reg = QRW_Reg(dim).cl
        self.search = state_search
        circ = self.circuit
        # Hadamards on all of the states
        circ.h(quantum_reg[0:self.n])
        # Operating U_prime pi/2 sqrt(2) times
        times = np.ceil(0.5 * np.pi * np.sqrt(2**(self.dim)))
        for i in range(int(times)):
            circ.unitary(self.U_prime(), quantum_reg, label="U'")
        # Measure the registry onto the classical bits
        circ.measure(quantum_reg, classical_reg)

    # Define all of the operators for the QRW search: S, C, U, U'

    def S(self, display_matrix=False):
        coin = self.c
        num_dir = int(2**coin)
        state = self.dim
        S = np.zeros((2**(coin + state), 2**(coin + state)))
        for i in range(num_dir):
          for j in range(2**state):
            # Perform the bit flip using a XOR
            j_bin = int(j)
            e_i = int(i + 1)
            xor = j_bin ^ e_i
            row = xor + (i * (2**state))
            S[row][j + (i * (2**state))] = 1
        if display_matrix:
          print('Matrix for the shift operator:')
          print(S)
        return Operator(S)

    def C(self, display_matrix=False):
        # Definition of the C operator as an outer product in the coin space --> tensor product with the identity in state space
        coin = self.c
        num_dir = int(2**coin)
        state = self.dim
        num_state = int(2**state)
        # Define the operator in just the coin space
        s_c = np.zeros((2**(coin), 2**(coin)))
        I_c = np.zeros((2**(coin), 2**(coin)))
        for i in range(num_dir):
            I_c[i][i] = 1
            for j in range(num_dir):
                s_c[i][j] = num_dir**(-1)
        s_c = 2 * s_c
        s_c = Operator(s_c)
        I_c = Operator(I_c)
        G = s_c - I_c
        # Define the identity in the state space
        I_s = np.zeros((2**(state), 2**(state)))
        for i in range(num_state):
            I_s[i][i] = 1
        I = Operator(I_s)
        C = G.tensor(I)
        if display_matrix:
            print('Matrix for the quantum coin operator:')
            print(np.real(C.data))
        return C

    def U(self, display_matrix=False):
        S_= self.S()
        C_ = self.C()
        U = C_.compose(S_)
        if display_matrix:
            print('Matrix for U:')
            print(np.real(U.data))
        return U

    def U_prime(self, display_matrix=False):
        # state_search is the state we are searching for --> focus on the second term
        # Note: the state search must be in decimal!
        if self.search >= 2**self.dim:
            raise TypeError("search_state parameter is outside of state space values.")
        elif self.search < 0:
            raise TypeError("search_state parameter is outside of state space values.")
        else:
            # Focus on the second term of U_prime
            coin = self.c
            num_dir = int(2**coin)
            state = self.dim
            num_state = int(2**state)
            search_array = np.zeros((num_state, num_state))
            search_array[self.search][self.search] = 1
            search = Operator(search_array)
            s_c = np.zeros((2**(coin), 2**(coin)))
            for i in range(num_dir):
                for j in range(num_dir):
                    s_c[i][j] = num_dir**(-1)
            coin_ = Operator(s_c)
            search_op = coin_.tensor(search)
            S_ = self.S()
            second_term = search_op.compose(S_)
            U_ = self.U()
            U_prime = U_ - (2 * second_term)
            if display_matrix:
                print("Matrix for U':")
                print(np.real(U_prime.data))
            return U_prime

    # Visualization
    def draw_circuit(self):
        return self.circuit.draw(output='mpl')

    # Plots by actually measuring the circuit
    def plot_states_hist(self):
        backend = Aer.get_backend('qasm_simulator')
        shots = 1024 # number of times circuit is run (for sampling)
        results = execute(self.circuit, backend=backend, shots=shots).result()
        counts = results.get_counts()
        return plot_histogram(counts, figsize=(5, 5), sort='value_desc')

# Define the QRW algorithm as a subclass of the QRW registry
class QRW(QRW_Reg):
    """
    The QRW class is an arbitrary implementation of a QRW
    dim: number of qubits in the state space
    c_label: string that determines what coin operator to use on the coin space
    - hadamard = Hadamard operator --> tensor product with the identity
    - random = random unitary operator
    step: number of 'steps' the QRW completes which corresponds to the number of times the operator U is called
    """
    def __init__(self, dim, c_label, step):
        QRW_Reg.__init__(self, dim)
        self.circuit = QRW_Reg(dim).circuit
        quantum_reg = QRW_Reg(dim).qr
        classical_reg = QRW_Reg(dim).cl
        circ = self.circuit
        # Hadamards on all of the states
        circ.h(quantum_reg[0:self.n])
        for i in range(int(step)):
            circ.unitary(self.U(c_label), quantum_reg, label="step")
        # Measure the registry onto the classical bits
        circ.measure(quantum_reg, classical_reg)

    # Defining all of the operators for the QRW search: S, C, U

    def S(self):
        coin = self.c
        num_dir = int(2**coin)
        state = self.dim
        S = np.zeros((2**(coin + state), 2**(coin + state)))
        for i in range(num_dir):
            for j in range(2**state):
                # Perform the bit flip using a XOR
                j_bin = int(j)
                e_i = int(i + 1)
                xor = j_bin ^ e_i
                row = xor + (i * (2**state))
                S[row][j + (i * (2**state))] = 1
        return Operator(S)

    def C(self, c_label):
        coin = self.c
        state = self.dim
        # Create the identity in the S space
        I = np.zeros((2**state, 2**state))
        for i in range(2**state):
            I[i][i] = 1
        I = Operator(I)
        if c_label == "hadamard":
            result= np.zeros((2**coin, 2**coin))
            for i in range(2**coin):
                for j in range(2**coin):
                    if i >= 2 and j >= 2:
                        result[i][j] = (-1) * (-1)**(i * j) * (2**((-1) * (0.5 * coin)))
                    else:
                        result[i][j] = (-1)**(i * j) * (2**((-1) * (0.5 * coin)))
            res_op = (Operator(result))
            C_final = res_op.tensor(I)
            return C_final
        elif c_label == "random":
            dim = []
            for i in range(coin):
                dim.append(2)
            res_op = random_unitary(tuple(dim))
            C_final = res_op.tensor(I)
            return C_final
        else:
            raise TypeError("Label string for C is not a valid input!")

    def U(self, c_label):
        S_= self.S()
        C_ = self.C(c_label)
        U = C_.compose(S_)
        return U

    # Visualization
    def draw_circuit(self):
        return self.circuit.draw()

    # Plot by actually measuring the circuit
    def plot_states_hist(self):
        backend = Aer.get_backend('qasm_simulator')
        print(backend)
        shots = 1024 # number of times circuit is run (for sampling)
        results = execute(self.circuit, backend=backend, shots=shots).result()
        counts = results.get_counts()
        return plot_histogram(counts, figsize=(5, 5), sort='value_desc')

    # Control the circuit
    def execute(self):
        backend = Aer.get_backend('qasm_simulator')
        results = execute(self.circuit, backend=backend, shots=1).result()
        answer = results.get_counts()
        # One execution means that there will be one state in the answer dictionary
        state = list(answer.keys())
        return state[0]

# Define the QRW search circuit
search_circuit = QRW_Search(2, 2)

search_circuit.draw_circuit()

search_circuit.plot_states_hist()

n_steps = 1
q = QuantumRegister(3, name='q') # 3 qubits register
c = ClassicalRegister(3, name='c') # 3 bits classical register
q_circ = QuantumCircuit(q, c) # quantum circuit
for i in range(n_steps):
    q_circ.h(0)
    q_circ.ccx(0, 2, 1)
    q_circ.ccx(1, 2, 0)
    q_circ.x(0)
    q_circ.ccx(0, 1, 2)
    q_circ.x(0)
    q_circ.x(2)
    q_circ.ccx(1, 2, 0)
    q_circ.ccx(0, 2, 1)
    q_circ.x(1)
    q_circ.ccx(1, 2, 0)
    q_circ.x(1)
    q_circ.x(2)
    q_circ.ccx(1, 2, 0)
    q_circ.x(0)
    q_circ.ccx(0, 1, 2)
    q_circ.x(0)
    q_circ.x(2)
    q_circ.ccx(1, 2, 0)
    q_circ.x(2)
    q_circ.ccx(0, 2, 1)
    q_circ.ccx(1, 2, 0)
    q_circ.x(0)
    q_circ.ccx(0, 1, 2)
    q_circ.x(0)
    q_circ.x(1)
    q_circ.ccx(1, 2, 0)
    q_circ.x(1)
    q_circ.ccx(0, 2, 1)
    q_circ.ccx(1, 2, 0)
q_circ.barrier()
q_circ.measure([q[0], q[1], q[2]], [c[0], c[1], c[2]])
q_circ.draw('mpl')

backend = Aer.get_backend('qasm_simulator')
shots = 1024 # number of times circuit is run (for sampling)
results = execute(q_circ, backend=backend, shots=shots).result()
counts = results.get_counts()
# Normalize counts to obtain probabilities
total_counts = sum(counts.values())
probabilities = {state: count / total_counts for state, count in counts.items()}
plot_histogram(probabilities, figsize=(5, 5), sort='value_desc')

n_steps = 4
q = QuantumRegister(3, name='q') # 3 qubits register
c = ClassicalRegister(3, name='c') # 3 bits classical register
q_circ = QuantumCircuit(q, c) # quantum circuit
for i in range(n_steps):
    q_circ.h(q[0]) # coin
    q_circ.ccx(q[0], q[1], q[2]) # Toffoli
    q_circ.cx (q[0], q[1]) # CNOT
    q_circ.x(q[0])
    q_circ.x(q[1])
    q_circ.ccx(q[0], q[1], q[2])
    q_circ.x(q[1])
    q_circ.cx(q[0], q[1])
q_circ.barrier()
q_circ.measure([q[0], q[1], q[2]], [c[0], c[1], c[2]])
q_circ.draw('mpl')

backend = Aer.get_backend('qasm_simulator')
shots = 1024 # number of times circuit is run (for sampling)
results = execute(q_circ, backend=backend, shots=shots).result()
counts = results.get_counts()
# Normalize counts to obtain probabilities
total_counts = sum(counts.values())
probabilities = {state: count / total_counts for state, count in counts.items()}
plot_histogram(probabilities, figsize=(5, 5), sort='value_desc')

# IBMQ.save_account needs to be run only once
IBMQ.save_account(API_KEY, overwrite=True)
# Load the saved IBMQ account
provider = IBMQ.load_account()
# Get a list of available backends
backend_overview()

# Get a specific backend device by its name
backend_device = provider.get_backend('ibm_brisbane')
print("Running on: "+str(backend_device))

# Monitor the status of the selected backend device
backend_monitor(backend_device)

shots = 1024 # number of times circuit is run (for sampling)
# Execute the quantum circuit on the selected backend
results_r = execute(q_circ, backend_device, shots=shots).result()
counts_r = results_r.get_counts()
# Normalize counts to obtain probabilities
total_counts = sum(counts_r.values())
probabilities_r = {state: count_r / total_counts for state, count_r in counts_r.items()}
plot_histogram([probabilities_r, probabilities], legend=['ibm_brisbane', 'qasm_simulator'], color=['#061727', '#82cfff'], sort='value_desc')

# Animate the 2D 'mapping' for the QRW it inherits the QRW class
class QRW_Automata(QRW):
    """
    QRW_Automata is a class that inherits the QRW class --> it animates multiple executions of the QRW algorithm into a cellular automata board
    Note: this algorithm has only been defined for the cases of 2 and 4 state qubits!
    dim: number of qubits in the state space
    c_label: string that determines what coin operator to use on the coin space
    - hadamard = Hadamard operator --> tensor product with the identity
    - random = random unitary operator
    steps: number of 'steps' the QRW completes which corresponds to the number of times the operator U is called
    iters: number of times the circuit is executed (also determines the number of frames in the animation)
    """
    def __init__(self, dim, c_label, steps, iters):
        QRW.__init__(self, dim, c_label, steps)
        self.n = QRW(dim, c_label, steps).n
        self.c = QRW(dim, c_label, steps).c
        self.circ = QRW(dim, c_label, steps).circuit
        state = []
        for i in range(iters):
            state.append(QRW(dim, c_label, steps).execute())
        self.state = state
        print("State:")
        print(state)
        c_state = []
        s_state = []
        for i in state:
            c_state.append(i[0:self.c])
            s_state.append(i[self.c:])
        self.c_state = c_state
        self.s_state = s_state
        if (dim/2).is_integer() != True:
            raise ValueError("The half of the number of qubits in the state space has to be an integer!")
        # Divide up the board according to number of bits
        # n: number of rows and columns --> we are going to implement only the case for n=4!
        if int(self.dim) != 4 and int(self.dim) != 2:
            raise ValueError("Grid is implemented only for the cases where n=2 and n=4 in the state space!")
        figure, axes = plt.subplots(nrows=2, ncols=1)
        axes[0].set_facecolor((0, 0, 0))
        axes[0].get_xaxis().set_ticks([])
        axes[0].get_yaxis().set_ticks([])
        axes[1].get_xaxis().set_ticks([])
        axes[1].get_yaxis().set_ticks([])
        self.fig = figure
        self.ax = axes
        self.circuit.draw('mpl', scale=0.9, ax=axes[1])
        n_v = (2 * self.dim)
        if self.c % 2 != 0:
          n_h = self.dim
        else:
          n_h = n_v
        v_lines = np.arange(0, 1, 1/n_v)
        h_lines = np.arange(0, 1, 1/n_h)
        # Draw the frames for the board
        for i in v_lines:
            axes[0].axvline(x=i, ymin=0, ymax=1, color='r')
        for i in h_lines:
            axes[0].axhline(y=i, xmin=0, xmax=1, color='r')
        anim = FuncAnimation(self.fig, self.animate, frames=int(len(state)), interval=1300)
        anim.save("QRW_{n}qubits_{op}_{st}steps.gif".format(n=self.dim, op=c_label, st=steps),)

    def animate(self, i):
        # Plot the state space
        if self.dim == 2:
            c_ = self.c_state[i]
            s_ = self.s_state[i]
            n_v = (2 * self.dim)
            n_h = self.dim
            verts = [
                (0. + int(s_[-1]) * (1/n_v) + 0.5 * int(c_), 1. - int(s_[0]) * (1/n_h) - (1/n_h)), # left, bottom
                (0. + int(s_[-1]) * (1/n_v)+ 0.5 * int(c_), 1. - int(s_[0]) * (1/n_h)), # left, top
                ((1/n_v) + int(s_[-1]) * (1/n_v) + 0.5 * int(c_), 1. - int(s_[0]) * (1/n_h)), # right, top
                ((1/n_v) + int(s_[-1]) * (1/n_v) + 0.5 * int(c_), 1. - int(s_[0]) * (1/n_h) - (1/n_h)), # right, bottom
                (0. + int(s_[-1]) * (1/n_v) + 0.5 * int(c_), 1. - int(s_[0]) * (1/n_h) - (1/n_h)), # ignored
            ]
            codes = [
                Path.MOVETO,
                Path.LINETO,
                Path.LINETO,
                Path.LINETO,
                Path.CLOSEPOLY,
            ]
            path = Path(verts, codes)
            cell = patch.PathPatch(path, facecolor='w')
            self.ax[0].set_xlabel("Measured state: {state_}".format(state_=self.state[i]))
            self.ax[0].add_patch(cell)
            return patch
        elif self.dim == 4:
            c_ = self.c_state[i]
            s_ = self.s_state[i]
            n_v = (2 * self.dim)
            n_h = (2 * self.dim)
            verts = [
                (0. + int(s_[-1]) * (1/n_v) + int(s_[-2]) * (2/n_v) + 0.5 * int(c_[1]), 1. - int(s_[0])* 2 *(1/n_h) - int(s_[1]) * (1/n_h) - (1/n_h) - 0.5 * int(c_[0])), # left, bottom
                (0. + int(s_[-1]) * (1/n_v) + int(s_[-2]) * (2/n_v) + 0.5 * int(c_[1]), 1. - int(s_[0])* 2 *(1/n_h) - int(s_[1]) * (1/n_h) - 0.5 * int(c_[0])), # left, top
                ((1/n_v) + int(s_[-1]) * (1/n_v) + int(s_[-2]) * (2/n_v) + 0.5 * int(c_[1]), 1. - int(s_[0])* 2 *(1/n_h) - int(s_[1]) * (1/n_h) - 0.5 * int(c_[0])), # right, top
                ((1/n_v) + int(s_[-1]) * (1/n_v) + int(s_[-2]) * (2/n_v) + 0.5 * int(c_[1]), 1. - int(s_[0])* 2 *(1/n_h) - int(s_[1]) * (1/n_h) - (1/n_h) - 0.5 * int(c_[0])), # right, bottom
                (0. + int(s_[-1]) * (1/n_v) + int(s_[-2]) * (2/n_v) + 0.5 * int(c_[1]), 1. - int(s_[0])* 2 *(1/n_h) - int(s_[1]) * (1/n_h) - (1/n_h) - 0.5 * int(c_[0])), # ignored
            ]
            codes = [
                Path.MOVETO,
                Path.LINETO,
                Path.LINETO,
                Path.LINETO,
                Path.CLOSEPOLY,
            ]
            path = Path(verts, codes)
            cell = patch.PathPatch(path, facecolor='w')
            self.ax[0].set_xlabel("Measured state: {state_}".format(state_=self.state[i]))
            self.ax[0].add_patch(cell)
            return patch

QRW_Automata(4, "hadamard", 4, 16)

QRW_Automata(4, "random", 4, 16)

display(Image(data=open('/content/QRW_4qubits_hadamard_4steps.gif','rb').read(), format='gif'))

display(Image(data=open('/content/QRW_4qubits_random_4steps.gif','rb').read(), format='gif'))