optimisation3_solution_quantum_challenge_2020.py

import typing as ty
from itertools import permutations

from qiskit import (Aer, ClassicalRegister, QuantumCircuit, QuantumRegister,
                    execute)
from qiskit.transpiler import PassManager
from qiskit.transpiler.passes import Unroller
from qprof import profile


def mcz(register_size: int) -> QuantumCircuit:
    """Returns a multi-controlled Z gate."""
    assert register_size > 1, "Cannot apply a multi-controlled Z without any control."
    ancilla_size = max(0, register_size - 2)
    qr = QuantumRegister(register_size)
    if ancilla_size > 0:
        qa = QuantumRegister(ancilla_size)
        circuit = QuantumCircuit(qr, qa, name="mcz")
        circuit.h(qr[-1])
        circuit.mcx(qr[:-1], qr[-1], qa, mode="v-chain")
        circuit.h(qr[-1])
    else:
        circuit = QuantumCircuit(qr, name="mcz")
        circuit.cz(qr[:-1], qr[-1])
    return circuit


def U_amp(qubit_numbers: int) -> QuantumCircuit:
    """Apply the amplification step of Grover's algorithm."""
    mcz_gate = mcz(qubit_numbers)

    q = QuantumRegister(qubit_numbers)
    if qubit_numbers > 2:
        qancill = QuantumRegister(qubit_numbers - 2)
        circuit = QuantumCircuit(q, qancill, name="U_amp")
    else:
        qancill = list()
        circuit = QuantumCircuit(q, name="U_amp")
    circuit.h(q)
    circuit.x(q)
    circuit.append(mcz_gate, [*q, *qancill])
    circuit.x(q)
    circuit.h(q)
    return circuit


BoardEncoding = ty.List[ty.Tuple[str, str]]


def encode_problem(problem: BoardEncoding) -> QuantumCircuit:
    """
    Encode the given problem on the qboard register if the
    qstate register is in the |1> state.

    This function returns a quantum circuit that encodes the given problem
    on the |qboard> register only if the |qstate> register is |1>.

    :param problem: a list of asteroid coordinates, given as tuples
        of strings.
    """
    qboard = QuantumRegister(16)
    qstate = QuantumRegister(1)
    circuit = QuantumCircuit(qboard, qstate, name="encode_problem")

    for x, y in problem:
        index = 4 * int(x) + int(y)
        circuit.cx(qstate[0], qboard[index])
    return circuit


def qram(problem_set: ty.List[BoardEncoding]):
    """
    Encodes the given problem set into the |qboard> register.

    The returned routine takes as inputs:
    - |qindex> (called |x> in the implementation) that encodes the indices for
      which we want to encode the boards. The |index> register will likely be
      given in a superposition of all the possible states.
    - |qboard> that should be given in the |0> state and that is returned in
      the state that encodes the problem_set given.
      For each qindex, the |qboard> state in |qindex>|qboard> is full of |0>
      except on the locations where there is an asteroid, in which case the
      state is |1>.
    - |qancilla> (called |q> in the implementation) is a 3-qubit ancilla
      register.

    This implementation is using  the "Select V" algorithm from
    https://arxiv.org/abs/1711.10980. See Section G.4. and Figure 6.

    Note: everything has been unrolled because I did not have the time to write
    and debug the generic version of the algorithm. Sorry for this.
    """

    x = QuantumRegister(4)  # Index
    qboard = QuantumRegister(16)
    q = QuantumRegister(3)
    circuit = QuantumCircuit(x, qboard, q, name="QRAM")

    # Trick to have the same indices as in the paper
    x = [None, *x[::-1]]
    q = [1, x[1], *q]
    # First, x[0] positive
    circuit.rccx(q[1], x[2], q[2])
    circuit.rccx(q[2], x[3], q[3])
    circuit.rccx(q[3], x[4], q[4])
    # q[4] is |1> only if the index was |1111>
    circuit.append(encode_problem(problem_set[0b1111]), [*qboard, q[4]])
    # Short-circuit (red)
    circuit.cx(q[3], q[4])
    # q[4] is |1> only if the index was |1110>
    circuit.append(encode_problem(problem_set[0b1110]), [*qboard, q[4]])
    # Short-circuit (green)
    circuit.cx(q[3], q[4])
    circuit.rccx(q[2], x[4], q[4])
    circuit.cx(q[2], q[3])
    # q[4] is |1> only if the index was |1101>
    circuit.append(encode_problem(problem_set[0b1101]), [*qboard, q[4]])
    # Short-circuit (red)
    circuit.cx(q[3], q[4])
    # q[4] is |1> only if the index was |1100>
    circuit.append(encode_problem(problem_set[0b1100]), [*qboard, q[4]])
    # Go up
    circuit.x(x[4])
    circuit.rccx(q[3], x[4], q[4])
    circuit.x(x[4]).inverse()
    # Short-circuit (green)
    circuit.cx(q[2], q[3])
    circuit.rccx(q[1], x[3], q[3])
    circuit.cx(q[1], q[2])
    # Go down
    circuit.rccx(q[3], x[4], q[4])
    # q[4] is |1> only if the index was |1011>
    circuit.append(encode_problem(problem_set[0b1011]), [*qboard, q[4]])
    # Short-circuit (red)
    circuit.cx(q[3], q[4])
    # q[4] is |1> only if the index was |1010>
    circuit.append(encode_problem(problem_set[0b1010]), [*qboard, q[4]])
    # Short-circuit (green)
    circuit.cx(q[3], q[4])
    circuit.rccx(q[2], x[4], q[4])
    circuit.cx(q[2], q[3])
    # q[4] is |1> only if the index was |1001>
    circuit.append(encode_problem(problem_set[0b1001]), [*qboard, q[4]])
    # Short-circuit (red)
    circuit.cx(q[3], q[4])
    # q[4] is |1> only if the index was |1000>
    circuit.append(encode_problem(problem_set[0b1000]), [*qboard, q[4]])
    # Go up
    circuit.x(x[4])
    circuit.rccx(q[3], x[4], q[4])
    circuit.x(x[4]).inverse()
    # Go up
    circuit.x(x[3])
    circuit.rccx(q[2], x[3], q[3])
    circuit.x(x[3]).inverse()
    # Short-circuit (green)
    # Small trick here: we are supposed to use q[0] but q[0] is "virtual"
    # and always equal to "1", so we can simplify.
    circuit.cx(q[1], q[2])
    circuit.cx(x[2], q[2])  # circuit.rccx(q[0], x[2], q[2])
    circuit.x(q[1])  # circuit.cx(q[0], q[1])
    # Go down
    circuit.rccx(q[2], x[3], q[3])
    # Go down
    circuit.rccx(q[3], x[4], q[4])
    # q[4] is |1> only if the index was |0111>
    circuit.append(encode_problem(problem_set[0b0111]), [*qboard, q[4]])
    # Short-circuit (red)
    circuit.cx(q[3], q[4])
    # q[4] is |1> only if the index was |0110>
    circuit.append(encode_problem(problem_set[0b0110]), [*qboard, q[4]])
    # Short-circuit (green)
    circuit.cx(q[3], q[4])
    circuit.rccx(q[2], x[4], q[4])
    circuit.cx(q[2], q[3])
    # q[4] is |1> only if the index was |0101>
    circuit.append(encode_problem(problem_set[0b0101]), [*qboard, q[4]])
    # Short-circuit (red)
    circuit.cx(q[3], q[4])
    # q[4] is |1> only if the index was |0100>
    circuit.append(encode_problem(problem_set[0b0100]), [*qboard, q[4]])
    # Go up
    circuit.x(x[4])
    circuit.rccx(q[3], x[4], q[4])
    circuit.x(x[4]).inverse()
    # Short-circuit (green)
    circuit.cx(q[2], q[3])
    circuit.rccx(q[1], x[3], q[3])
    circuit.cx(q[1], q[2])
    # Go down
    circuit.rccx(q[3], x[4], q[4])
    # q[4] is |1> only if the index was |0011>
    circuit.append(encode_problem(problem_set[0b0011]), [*qboard, q[4]])
    # Short-circuit (red)
    circuit.cx(q[3], q[4])
    # q[4] is |1> only if the index was |0010>
    circuit.append(encode_problem(problem_set[0b0010]), [*qboard, q[4]])
    # Short-circuit (green)
    circuit.cx(q[3], q[4])
    circuit.rccx(q[2], x[4], q[4])
    circuit.cx(q[2], q[3])
    # q[4] is |1> only if the index was |0001>
    circuit.append(encode_problem(problem_set[0b0001]), [*qboard, q[4]])
    # Short-circuit (red)
    circuit.cx(q[3], q[4])
    # q[4] is |1> only if the index was |0001>
    circuit.append(encode_problem(problem_set[0b0000]), [*qboard, q[4]])
    # Go up
    circuit.x(x[4])
    circuit.rccx(q[3], x[4], q[4])
    circuit.x(x[4]).inverse()
    # Go up
    circuit.x(x[3])
    circuit.rccx(q[2], x[3], q[3])
    circuit.x(x[3]).inverse()
    # Go up
    circuit.x(x[2])
    circuit.rccx(q[1], x[2], q[2])
    circuit.x(x[2]).inverse()
    # To avoid changing x[1]:
    circuit.x(x[1])
    return circuit


def check_if_match_permutation(permutation: ty.List[int]) -> QuantumCircuit:
    """
    Construct a quantum circuit to check if the given qboard state
    "includes" the given permutation.

    The |qresult> qubit is flipped if |qboard> includes the given
    permutation.
    """
    qboard = QuantumRegister(16)
    qresult = QuantumRegister(1)
    qancilla = QuantumRegister(2)
    circuit = QuantumCircuit(
        qboard, qresult, qancilla, name="check_if_match_permutation"
    )

    qubit_indices = [4 * row + col for row, col in enumerate(permutation)]
    circuit.mcx(
        [qboard[i] for i in qubit_indices],
        qresult,
        ancilla_qubits=qancilla,
        mode="v-chain",
    )

    return circuit


def get_2_qubit_adder_gate():
    """
    Returns a quantum circuit implementing an adder.

    The adder adds a 1-qubit value to a 2-bit register without
    any overflow check.
    """
    a = QuantumRegister(1)
    b = QuantumRegister(2)
    circuit = QuantumCircuit(a, b, name="2_qubit_adder_gate")

    circuit.ccx(a, b[0], b[1])
    circuit.cx(a, b[0])

    return circuit


def check_if_match_any_permutation():
    """
    Flip the phase of the boards that include at least 1 permutation.

    This quantum circuit takes a board on 16 qubits and flip the phase
    of the quantum state if the board includes at least 1 permutation.
    """
    qboard = QuantumRegister(16)
    qadded = QuantumRegister(1)
    qcounter = QuantumRegister(2)
    qancilla = QuantumRegister(2)
    circuit = QuantumCircuit(
        qboard, qadded, qcounter, qancilla, name="check_if_match_any_permutation"
    )

    add_gate = get_2_qubit_adder_gate()
    all_permutations = list(permutations(range(4), 4))

    for perm in all_permutations:
        checker = check_if_match_permutation(perm)
        circuit.append(checker, [*qboard, *qadded, *qancilla])
        circuit.append(add_gate, [*qadded, *qcounter])
        circuit.append(checker.inverse(), [*qboard, *qadded, *qancilla])

    circuit.z(qcounter)

    for perm in reversed(all_permutations):
        checker = check_if_match_permutation(perm)
        circuit.append(checker.inverse(), [*qboard, *qadded, *qancilla])
        circuit.append(add_gate.inverse(), [*qadded, *qcounter])
        circuit.append(checker, [*qboard, *qadded, *qancilla])

    return circuit


def grover_oracle_week3(problem_set) -> QuantumCircuit:
    """
    Flip the phase of the only index representing a non-solvable board.

    This oracle needs a total of 23 qubits:
    - 4 qubits for the |index> register
 qubits for the |board> register, that should be given in the
      |0> state.
    - 3 ancilla qubits.
    """
    qindex = QuantumRegister(4)
    qboard = QuantumRegister(16)
    qancilla = QuantumRegister(5)
    circuit = QuantumCircuit(qindex, qboard, qancilla, name="grover_oracle_week3")

    qram_gate = qram(problem_set)
    checker = check_if_match_any_permutation()

    circuit.append(qram_gate, [*qindex, *qboard, *qancilla[:3]])
    circuit.append(checker, [*qboard, *qancilla])
    circuit.append(qram_gate.inverse(), [*qindex, *qboard, *qancilla[:3]])

    return circuit


def week3_ans_func(problem_set) -> QuantumCircuit:
    """
    Returns a quantum circuit solving the given problem.
    """
    qindex = QuantumRegister(4)
    qboard = QuantumRegister(16)
    qancilla = QuantumRegister(5)
    cindex = ClassicalRegister(4)
    circuit = QuantumCircuit(qindex, qboard, qancilla, cindex, name="week3_ans_func")

    oracle = grover_oracle_week3(problem_set)
    diffusion = U_amp(4)

    circuit.h(qindex)

    # Perform only 1 iteration to stick to the rules.
    for _ in range(1):
        circuit.append(oracle, [*qindex, *qboard, *qancilla])
        circuit.append(diffusion, [*qindex, *qancilla[:2]])

    circuit.measure(qindex, cindex)
    circuit = circuit.reverse_bits()
    return circuit


def cost(circuit: QuantumCircuit) -> int:
    """Compute the cost of a circuit"""
    pass_ = Unroller(["u3", "cx"])
    pm_unroll = PassManager(pass_)
    unrolled_circuit = pm_unroll.run(circuit)
    ops = unrolled_circuit.count_ops()
    return 10 * ops.get("cx", 0) + ops.get("u3", 0)


problem_set = [
    [["0", "2"], ["1", "0"], ["1", "2"], ["1", "3"], ["2", "0"], ["3", "3"]],
    [["0", "0"], ["0", "1"], ["1", "2"], ["2", "2"], ["3", "0"], ["3", "3"]],
    [["0", "0"], ["1", "1"], ["1", "3"], ["2", "0"], ["3", "2"], ["3", "3"]],
    [["0", "0"], ["0", "1"], ["1", "1"], ["1", "3"], ["3", "2"], ["3", "3"]],
    [["0", "2"], ["1", "0"], ["1", "3"], ["2", "0"], ["3", "2"], ["3", "3"]],
    [["1", "1"], ["1", "2"], ["2", "0"], ["2", "1"], ["3", "1"], ["3", "3"]],
    [["0", "2"], ["0", "3"], ["1", "2"], ["2", "0"], ["2", "1"], ["3", "3"]],
    [["0", "0"], ["0", "3"], ["1", "2"], ["2", "2"], ["2", "3"], ["3", "0"]],
    [["0", "3"], ["1", "1"], ["1", "2"], ["2", "0"], ["2", "1"], ["3", "3"]],
    [["0", "0"], ["0", "1"], ["1", "3"], ["2", "1"], ["2", "3"], ["3", "0"]],
    [["0", "1"], ["0", "3"], ["1", "2"], ["1", "3"], ["2", "0"], ["3", "2"]],
    [["0", "0"], ["1", "3"], ["2", "0"], ["2", "1"], ["2", "3"], ["3", "1"]],
    [["0", "1"], ["0", "2"], ["1", "0"], ["1", "2"], ["2", "2"], ["2", "3"]],
    [["0", "3"], ["1", "0"], ["1", "3"], ["2", "1"], ["2", "2"], ["3", "0"]],
    [["0", "2"], ["0", "3"], ["1", "2"], ["2", "3"], ["3", "0"], ["3", "1"]],
    [["0", "1"], ["1", "0"], ["1", "2"], ["2", "2"], ["3", "0"], ["3", "1"]],
]

circuit = week3_ans_func(problem_set)

print(f"Cost of the generated quantum circuit: {cost(circuit)}")