qiskit exercise on Grover search algorithm

exercises/grover_exercise.py

+#!/bin/env python3
+
+'''
+Exercise with Qiskit: Grover algorithm.
+
+Task: Implement the three missing parts marked by "TODO" in the function grover().
+
+This script requires qiskit. Further information about Qiskit can be found
+here: https://qiskit.org
+
+Installation instructions for Qiskit and links to cloud environments
+with Qiskit pre-installed can be found here:
+https://qiskit.org/documentation/getting_started.html
+
+Running this script does not require a full installation of Qiskit.
+Installing qiskit-terra should be sufficient:
+https://github.com/Qiskit/qiskit-terra
+However, with qiskit-aer also installed the simulation will be much faster:
+https://github.com/Qiskit/qiskit-aer
+'''
+
+###############################################################################
+###   LOAD MODULES
+###############################################################################
+
+from traceback import print_exc
+from sys import stderr
+import numpy as np # numpy is a dependency of qiskit and should be installed.
+
+# Import basic functions and classes of Qiskit.
+# The functions/classes are imported explicitly to obtain more useful error
+# messages if this fails.
+try:
+    from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, execute
+    from qiskit.circuit.library import XGate
+except:
+    print("Failed to load qiskit:", file=stderr)
+    print_exc(file=stderr)
+    print(__doc__)
+    exit(1)
+
+# Load a simulation backend. Depending on the Qiskit installation, either
+# a simulator from qiskit-aer or a basic simulator from qiskit-terra is used.
+try:
+    # If qiskit-aer is installed: use the backends from this package.
+    from qiskit import Aer
+    backend = Aer.get_backend('qasm_simulator')
+    MAX_QUBITS = 14
+except:
+    # If qiskit-aer is can not be loaded: use the basic backends provided by qiskit-terra.
+    from qiskit import BasicAer
+    backend = BasicAer.get_backend('qasm_simulator')
+    # The basic backend is slow compared to the backend in Aer.
+    print('WARNING: using the slow simulation backend from qiskit-terra.\nFaster simulations with more qubits and shots are possible if qiskit-aer is installed.')
+    MAX_QUBITS = 9
+
+### SIMULATE REAL DEVICE
+# If you want to simulate a real quantum device, you can use this backend:
+# <https://qiskit.org/documentation/tutorials/simulators/2_device_noise_simulation.html>
+#from qiskit.test.mock import FakeVigo
+#backend = FakeVigo()
+#MAX_QUBITS = 5
+
+###############################################################################
+
+
+def random_bitstring(n):
+    '''
+    Create a random string of length n containing only 0s and 1s.
+    '''
+    bitstr = bin(np.random.randint(2**n))[2:]
+    bitstr = (n-len(bitstr))*'0' + bitstr
+    return bitstr
+
+
+def grover(
+        n : 'number of bits' = 3,
+        shots : 'number of simulation runs' = 1000,
+        show_circuit : 'print quantum circuits to stdout' = True,
+        ):
+    '''
+    Simulate the Grover search algorithm.
+    '''
+    # Don't try to use too many qubits.
+    if n >= MAX_QUBITS:
+        raise ValueError("Circuit is too large")
+
+    # Create a random bit string <solution>.
+    solution = random_bitstring(n)
+    # The following string may be useful later on.
+    zeros = n*'0'
+
+    print('Generated random bit string: %s'%solution)
+
+
+    ### INITIALIZATION ###
+
+    # Create a register of n qubits.
+    qubit_register = QuantumRegister(n, name='q')
+    # Create one ancilla qubit.
+    ancilla = QuantumRegister(1, name='a')
+    # Create a circuit containing the previously created registers.
+    full_circuit = QuantumCircuit(qubit_register, ancilla)
+
+    # All qubits are now initialized in the |0> state.
+
+    # Apply a Hadamard gate on each of the qubits in <qubit_register>.
+    full_circuit.h(qubit_register)
+
+    # TODO: bring the ancilla qubit in the |-> state.
+    # HINT: You can use Pauli gates (full_circuit.x, full_circuit.y,
+    #       full_circuit.z) and the Hadamard gate (full_circuit.h).
+
+    if show_circuit:
+        print('Initialization before the Grover steps:')
+        print(full_circuit.draw('text'))
+
+
+    ### GROVER STEP ###
+
+    # Create a new quantum circuit with the same qubits, in which we implement
+    # a single Grover step. We will later include this circuit multiple times
+    # in <full_circuit>.
+    step_circuit = QuantumCircuit(qubit_register, ancilla, name='G')
+
+    # The oracle is defined as follows:
+    # If and only if the qubits in <qubit_register> are in the state
+    # <solution>, a Pauli X gate will be applied to the ancilla qubit.
+    # The following command creates such a gate (<n> is the number of
+    # control qubits):
+    oracle_gate = XGate().control(n, ctrl_state=solution)
+
+    # Apply the oracle gate to the Grover step circuit.
+    # Here the second argument specifies on which qubits the gate acts.
+    # In this case, these are all qubits (qubit_register and ancilla).
+    all_qubits = (*qubit_register, *ancilla)
+    step_circuit.append(oracle_gate, all_qubits)
+
+    # TODO: Implement the remaining part of the Grover step.
+    # HINT: You can ignore global phases. The previously introduced functions
+    #       and variables are sufficient to solve the task.
+
+    # Show the quantum circuit for a single Grover step.
+    if show_circuit:
+        print('Circuit for a single Grover step G:')
+        print(step_circuit.draw('text'))
+
+    # TODO: Calculate how many Grover steps (iterations) should be applied and
+    #       correct the following line.
+    # HINT: You can use round(...) to obtain an integer.
+    #       "2 to the power n" is implemented as 2**n in python.
+    iterations = 1
+    # Add the grover step circuit <iteration> times to the full circuit.
+    full_circuit.compose( step_circuit.power(iterations), inplace=True )
+
+
+    ### MEASUREMENT ###
+
+    # Create a classical register of n bits.
+    cbit_register = ClassicalRegister(n, 'c')
+    # Add <cbit_register> to the quantum circuit.
+    full_circuit.add_register(cbit_register)
+    # Measure the qubits in <qubit_register> and save the values to
+    # <cbit_register>.
+    full_circuit.measure(qubit_register, cbit_register)
+
+
+    # Show the (now complete) quantum circuit.
+    if show_circuit:
+        print('Full quantum circuit:')
+        print(full_circuit.draw('text'))
+
+
+    ### SIMULATION ###
+
+    # Create a simulation job for this circuit.
+    job = execute(full_circuit, backend=backend, shots=shots)
+    result = job.result()
+    # Get the outputs of the simulation.
+    # counts is a dictionary mapping the classical outputs to the number
+    # how often this output occurred. The classical output is given as a
+    # bit string with bits ordered from bit n to 0.
+    counts = result.get_counts()
+
+
+    ### SHOW RESULTS ###
+
+    # Show the output of the simulation.
+    for key, value in counts.items():
+        print('Measured %s in%7d of %d simulation runs.'%(key, value, shots))
+    # Show the input bit string  a  for comparison.
+    print('Check:   %s is the correct solution.'%(''.join(str(i) for i in solution)))
+
+    # Check if the circuit solved the problem of finding the right solution
+    # with high probability.
+    if counts[solution] > 0.9 * shots:
+        print('It works: correct result obtained in {}% of the simulation runs.'\
+                .format(100*counts[solution]/shots))
+    elif counts[solution] < 0.6 * shots:
+        print('Improvement required: correct result obtained in only {}% of the simulation runs.'\
+                .format(100*counts[solution]/shots))
+        print('Implement the three parts marked by "TODO" in the function grover().')
+
+    # Analyze the success probability.
+    # Instead of p_success we print p_fail, because for large n and small p_fail
+    # that gives a better overview.
+    p_fail_theory = np.cos((2*iterations + 1)*np.arccos((1-2**-n)**.5))**2
+    p_fail_simulation = (shots - counts[solution])/shots
+    print('Probability of failure: theory = %g,  simulation = %g'\
+            %(p_fail_theory, p_fail_simulation))
+
+
+
+# This is how a python script should be structured:
+# In the end the main function is called if the script is executed directly.
+if __name__ == '__main__':
+    grover()

exercises/grover_mpl.py

+#!/usr/bin/env python3
+
+'''
+Exercise with Qiskit: Grover algorithm.
+
+Task: Implement the three missing parts marked by "TODO" in the function grover().
+
+This script requires qiskit. Further information about Qiskit can be found
+here: https://qiskit.org
+
+Installation instructions for Qiskit and links to cloud environments
+with Qiskit pre-installed can be found here:
+https://qiskit.org/documentation/getting_started.html
+
+Running this script does not require a full installation of Qiskit.
+Installing qiskit-terra should be sufficient:
+https://github.com/Qiskit/qiskit-terra
+However, with qiskit-aer also installed the simulation will be much faster:
+https://github.com/Qiskit/qiskit-aer
+'''
+
+###############################################################################
+###   LOAD MODULES
+###############################################################################
+
+from traceback import print_exc
+from sys import stderr
+import numpy as np
+import matplotlib.pyplot as plt
+
+from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, execute
+from qiskit.circuit.library import XGate, ZGate
+
+# Load a simulation backend. Depending on the Qiskit installation, either
+# a simulator from qiskit-aer or a basic simulator from qiskit-terra is used.
+try:
+    # If qiskit-aer is installed: use the backends from this package.
+    from qiskit import Aer
+    default_backend = Aer.get_backend('qasm_simulator')
+    MAX_QUBITS = 14
+except:
+    # If qiskit-aer is can not be loaded: use the basic backends provided by qiskit-terra.
+    from qiskit import BasicAer
+    default_backend = BasicAer.get_backend('qasm_simulator')
+    # The basic backend is slow compared to the backend in Aer.
+    print('WARNING: using the slow simulation backend from qiskit-terra.\nFaster simulations with more qubits and shots are possible if qiskit-aer is installed.')
+    MAX_QUBITS = 9
+
+###############################################################################
+
+def to_bit_string(x, n):
+    bitstr = bin(x)[2:]
+    return (n-len(bitstr))*'0' + bitstr
+
+def random_bitstring(n):
+    '''
+    Create a random string of length n containing only 0s and 1s.
+    '''
+    return to_bit_string(np.random.randint(2**n), n)
+
+def visualize(circuit, title, fmt):
+    try:
+        visualization = circuit.draw(fmt)
+    except:
+        print('Failed to visualize circuit')
+        return
+    if fmt == 'mpl':
+        visualization.text(0.5, 0.9, title, horizontalalignment='center')
+    else:
+        print(title)
+        print(visualization)
+
+
+def grover(
+        n : 'number of bits' = 3,
+        shots : 'number of simulation runs' = 1000,
+        iterations : 'number of Grover steps or None' = None,
+        show_circuit : 'format for showng quantum circuits' = 'text',
+        solution = None,
+        backend = default_backend,
+        ):
+    '''
+    Simulate the Grover search algorithm.
+    '''
+    # Don't try to use too many qubits.
+    if n >= MAX_QUBITS:
+        raise ValueError("Circuit is too large")
+
+    # Create a random bit string <solution>.
+    if solution is None:
+        solution = random_bitstring(n)
+    elif type(solution) == int:
+        solution = to_bit_string(solution, n)
+    # The following string may be useful later on.
+    zeros = n*'0'
+
+    print('Generated random bit string: %s'%solution)
+
+
+    ### INITIALIZATION ###
+
+    # Create a register of n qubits.
+    qubit_register = QuantumRegister(n, name='q')
+    # Create one ancilla qubit.
+    ancilla = QuantumRegister(1, name='a')
+    # Create a circuit containing the previously created registers.
+    full_circuit = QuantumCircuit(qubit_register, ancilla)
+
+    # All qubits are now initialized in the |0> state.
+
+    # Apply a Hadamard gate on each of the qubits in <qubit_register>.
+    full_circuit.h(qubit_register)
+
+    # TODO: bring the ancilla qubit in the |-> state.
+    # BEGIN SOLUTION
+    full_circuit.x(ancilla) # X|0> = |1>
+    full_circuit.h(ancilla) # H|1> = |->
+    # Alternative solution:
+    #full_circuit.h(ancilla) # H|0> = |+>
+    #full_circuit.z(ancilla) # Z|+> = |->
+    # END SOLUTION
+
+    if show_circuit:
+        visualize(full_circuit, 'Initialization before the Grover steps', show_circuit)
+
+
+    ### GROVER STEP ###
+
+    # Create a new quantum circuit with the same qubits, in which we implement
+    # a single Grover step. We will later include this circuit multiple times
+    # in <full_circuit>.
+    step_circuit = QuantumCircuit(qubit_register, ancilla, name='G')
+
+    # The oracle is defined as follows:
+    # If and only if the qubits in <qubit_register> are in the state
+    # <solution>, a Pauli X gate will be applied to the ancilla qubit.
+    # The following command creates such a gate (<n> is the number of
+    # control qubits):
+    oracle_gate = XGate().control(n, ctrl_state=solution)
+
+    # Apply the oracle gate to the Grover step circuit.
+    # Here the second argument specifies on which qubits the gate acts.
+    # In this case, these are all qubits (qubit_register and ancilla).
+    all_qubits = (*qubit_register, *ancilla)
+    step_circuit.append(oracle_gate, all_qubits)
+
+    # TODO: Implement the remaining part of the Grover step.
+    # HINT: You can ignore global phases. The previously introduced functions
+    #       and variables are sufficient to solve the task.
+    # BEGIN SOLUTION
+    step_circuit.h(qubit_register[:-1])
+    step_circuit.z(qubit_register[-1])
+    step_circuit.append(XGate().control(n-1, ctrl_state=(n-1)*'0'), qubit_register)
+    step_circuit.z(qubit_register[-1])
+    step_circuit.h(qubit_register[:-1])
+    # Alternative solution:
+    #step_circuit.h(qubit_register)
+    #step_circuit.barrier()
+    #step_circuit.x(qubit_register)
+    #step_circuit.mcrz(np.pi, qubit_register[1:], qubit_register[0])
+    #step_circuit.x(qubit_register)
+    #step_circuit.barrier()
+    #step_circuit.h(qubit_register)
+    # END SOLUTION
+
+    # Show the quantum circuit for a single Grover step.
+    if show_circuit:
+        visualize(step_circuit, 'Circuit for a single Grover step G', show_circuit)
+
+    # TODO: Calculate how many Grover steps (iterations) should be applied and
+    #       correct the following line.
+    # BEGIN SOLUTION
+    if (iterations is None):
+        #iterations = int(np.pi/(2*np.arcsin(2**(-n/2))))
+        iterations = int(np.pi*2**(n/2-2))
+    # END SOLUTION
+    # Add the grover step circuit <iteration> times to the full circuit.
+    full_circuit.append(step_circuit.power(iterations), all_qubits)
+
+
+    ### MEASUREMENT ###
+
+    # Create a classical register of n bits.
+    cbit_register = ClassicalRegister(n, 'c')
+    # Add <cbit_register> to the quantum circuit.
+    full_circuit.add_register(cbit_register)
+    # Measure the qubits in <qubit_register> and save the values to
+    # <cbit_register>.
+    full_circuit.measure(qubit_register, cbit_register)
+
+
+    # Show the (now complete) quantum circuit.
+    if show_circuit:
+        visualize(full_circuit, 'Full quantum circuit', show_circuit)
+
+
+    ### SIMULATION ###
+
+    # Create a simulation job for this circuit.
+    job = execute(full_circuit, backend=backend, shots=shots)
+    result = job.result()
+    # Get the outputs of the simulation.
+    # counts is a dictionary mapping the classical outputs to the number
+    # how often this output occurred. The classical output is given as a
+    # bit string with bits ordered from bit n to 0.
+    counts = result.get_counts()
+
+
+    ### SHOW RESULTS ###
+
+    # Show the output of the simulation.
+    for key, value in counts.items():
+        print('Measured %s in%7d of %d simulation runs.'%(key, value, shots))
+    # Show the input bit string  a  for comparison.
+    print('Check:   %s is the correct solution.'%(''.join(str(i) for i in solution)))
+
+    # Check if the circuit solved the problem of finding the right solution
+    # with high probability.
+    if counts[solution] > 0.9 * shots:
+        print('It works: correct result obtained in {}% of the simulation runs.'\
+                .format(100*counts[solution]/shots))
+    elif counts[solution] < 0.6 * shots:
+        print('Improvement required: correct result obtained in only {}% of the simulation runs.'\
+                .format(100*counts[solution]/shots))
+        print('Implement the three parts marked by "TODO" in the function grover().')
+
+    # Analyze the success probability.
+    # Instead of p_success we print p_fail, because for large n and small p_fail
+    # that gives a better overview.
+    p_fail_theory = np.cos((2*iterations + 1)*np.arcsin(2**(-n/2)))**2
+    p_fail_simulation = (shots - counts[solution])/shots
+    print('Probability of failure: theory = %g,  simulation = %g'\
+            %(p_fail_theory, p_fail_simulation))
+
+    return counts, solution
+
+
+def groverGUI(n=3, shots=1000, noisy=False, **kwargs):
+    if noisy:
+        from qiskit.test.mock import FakeVigo
+        backend = FakeVigo()
+    else:
+        backend = default_backend
+    counts, solution = grover(n=n, shots=shots, show_circuit='mpl', backend=backend, **kwargs)
+    solution = int(solution, 2)
+    fig, ax = plt.subplots()
+    ax.bar(range(2**n), [counts.get(to_bit_string(x, n), 0)/shots for x in range(2**n)], log=True)
+    ax.set_title('Simulation results')
+    ax.vlines(solution, 0, 1, color='red')
+    plt.show()
+
+
+# This is how a python script should be structured:
+# In the end the main function is called if the script is executed directly.
+if __name__ == '__main__':
+    from sys import argv
+    args = {}
+    for text in argv[1:]:
+        try:
+            key, value = text.split('=')
+            args[key] = int(value)
+        except:
+            print('Argument not understood:', text)
+    groverGUI(**args)