Neophytes team

.gitignore

@@ -1,4 +1,6 @@
 .venv
 .DS_Store
 .ipynb_checkpoints
-__pycache__
+__pycache__
+.idea
+dev

dev/maxcut.py

@@ -0,0 +1,176 @@
+import time
+from typing import Callable, List, Tuple
+
+import networkx as nx
+import numpy as np
+from pytket import Circuit, Qubit
+from pytket.backends.backend import Backend
+from pytket.backends.backendresult import BackendResult
+from pytket.extensions.qiskit import AerBackend
+from pytket.extensions.quantinuum import QuantinuumBackend
+from pytket.passes import DecomposeBoxes
+from pytket.pauli import Pauli, QubitPauliString
+from pytket.utils import QubitPauliOperator, gen_term_sequence_circuit
+
+from maxcut_plotting import plot_maxcut_results
+
+
+def qaoa_graph_to_cost_hamiltonian(edges: List[Tuple[int, int]], cost_angle: float) -> QubitPauliOperator:
+    qpo_dict = {QubitPauliString(): len(edges) * 0.5 * cost_angle}
+    for e in edges:
+        term_string = QubitPauliString([Qubit(e[0]), Qubit(e[1])], [Pauli.Z, Pauli.Z])
+        qpo_dict[term_string] = -0.5 * cost_angle
+    return QubitPauliOperator(qpo_dict)
+
+
+def qaoa_initial_circuit(n_qubits: int) -> Circuit:
+    c = Circuit(n_qubits)
+    for i in range(n_qubits):
+        c.H(i)
+    return c
+
+
+def qaoa_max_cut_circuit(
+    edges: List[Tuple[int, int]],
+    n_nodes: int,
+    mixer_angles: List[float],
+    cost_angles: List[float],
+) -> Circuit:
+    assert len(mixer_angles) == len(cost_angles)
+
+    # initial state
+    qaoa_circuit = qaoa_initial_circuit(n_nodes)
+
+    # add cost and mixer terms to state
+    for cost, mixer in zip(cost_angles, mixer_angles):
+        cost_ham = qaoa_graph_to_cost_hamiltonian(edges, cost)
+        mixer_ham = QubitPauliOperator({QubitPauliString([Qubit(i)], [Pauli.X]): mixer for i in range(n_nodes)})
+        qaoa_circuit.append(gen_term_sequence_circuit(cost_ham, Circuit(n_nodes)))
+        qaoa_circuit.append(gen_term_sequence_circuit(mixer_ham, Circuit(n_nodes)))
+
+    DecomposeBoxes().apply(qaoa_circuit)
+    return qaoa_circuit
+
+
+def max_cut_energy(edges: List[Tuple[int, int]], results: BackendResult) -> float:
+    energy = 0.0
+    dist = results.get_distribution()
+    for i, j in edges:
+        energy += sum((meas[i] ^ meas[j]) * prob for meas, prob in dist.items())
+
+    return energy
+
+
+def qaoa_instance_simple(
+    backend: Backend,
+    compiler_pass: Callable[[Circuit], bool],
+    guess_mixer_angles: np.array,
+    guess_cost_angles: np.array,
+    seed: int,
+    shots: int = 5000,
+) -> float:
+    # step 1: get state guess
+    my_prep_circuit = qaoa_max_cut_circuit(max_cut_graph_edges, n_nodes, guess_mixer_angles, guess_cost_angles)
+    measured_circ = my_prep_circuit.copy().measure_all()
+    compiler_pass(measured_circ)
+    res = backend.run_circuit(measured_circ, shots, seed=seed)
+
+    return max_cut_energy(max_cut_graph_edges, res)
+
+
+def qaoa_optimise_energy(
+    compiler_pass: Callable[[Circuit], bool],
+    backend: Backend,
+    iterations: int = 100,
+    n: int = 3,
+    shots: int = 5000,
+    seed: int = 12345,
+):
+    highest_energy = 0
+    best_guess_mixer_angles = [0 for i in range(n)]
+    best_guess_cost_angles = [0 for i in range(n)]
+    rng = np.random.default_rng(seed)
+    # guess some angles (iterations)-times and try if they are better than the best angles found before
+
+    for i in range(iterations):
+
+        guess_mixer_angles = rng.uniform(0, 1, n)
+        guess_cost_angles = rng.uniform(0, 1, n)
+
+        qaoa_energy = qaoa_instance_simple(
+            backend,
+            compiler_pass,
+            guess_mixer_angles,
+            guess_cost_angles,
+            seed=seed,
+            shots=shots,
+        )
+
+        if qaoa_energy > highest_energy:
+            print("new highest energy found: ", qaoa_energy)
+
+            best_guess_mixer_angles = guess_mixer_angles
+            best_guess_cost_angles = guess_cost_angles
+            highest_energy = qaoa_energy
+
+    print("highest energy: ", highest_energy)
+    print("best guess mixer angles: ", best_guess_mixer_angles)
+    print("best guess cost angles: ", best_guess_cost_angles)
+    return best_guess_mixer_angles, best_guess_cost_angles
+
+
+def qaoa_calculate(
+    backend: Backend,
+    compiler_pass: Callable[[Circuit], bool],
+    shots: int = 5000,
+    iterations: int = 100,
+    seed: int = 12345,
+) -> float:
+    # find the parameters for the highest energy
+    best_mixer, best_cost = qaoa_optimise_energy(compiler_pass, backend, iterations, 3, shots=shots, seed=seed)
+
+    # get the circuit with the final parameters of the optimisation:
+    my_qaoa_circuit = qaoa_max_cut_circuit(max_cut_graph_edges, n_nodes, best_mixer, best_cost)
+
+    my_qaoa_circuit.measure_all()
+
+    compiler_pass(my_qaoa_circuit)
+    handle = backend.process_circuit(my_qaoa_circuit, shots, seed=seed)
+
+    result = backend.get_result(handle)
+
+    return result
+
+
+max_cut_graph_edges = [(0, 1), (1, 2), (1, 3), (3, 4), (4, 5), (4, 6)]
+n_nodes = 7
+
+max_cut_graph = nx.Graph()
+max_cut_graph.add_edges_from(max_cut_graph_edges)
+# nx.draw(max_cut_graph, labels={node: node for node in max_cut_graph.nodes()})
+
+expected_results = [(0, 1, 0, 0, 1, 0, 0), (1, 0, 1, 1, 0, 1, 1)]
+
+cost_angle = 1.0
+cost_ham_qpo = qaoa_graph_to_cost_hamiltonian(max_cut_graph_edges, cost_angle)
+print(cost_ham_qpo)
+
+backend = AerBackend()
+# backend = QuantinuumBackend("H1-2E")
+# Total time for 100 iterations (ms): 87038.1588935852 80150.33507347107
+comp = backend.get_compiled_circuit
+iters = 10
+# 86
+start = time.time()
+res = qaoa_calculate(
+    backend,
+    backend.default_compilation_pass(2).apply,
+    shots=5000,
+    iterations=iters,
+    seed=12345,
+)
+
+end = time.time()
+print(f"Total time for {iters} iterations (ms): {(end - start) * 1000}")
+
+# plot_maxcut_results(res, 6)

dev/maxcut_plotting.py

@@ -0,0 +1,31 @@
+from pytket.backends.backendresult import BackendResult
+import matplotlib.pyplot as plt
+
+
+def plot_maxcut_results(result: BackendResult, n_strings: int) -> None:
+    """
+    Plots Maxcut results in a barchart with the two most common bitstrings highlighted in green.
+    """
+    counts_dict = result.get_counts()
+    sorted_shots = counts_dict.most_common()
+    n_shots = sum(counts_dict.values())
+
+    n_most_common_strings = sorted_shots[:n_strings]
+    x_axis_values = [str(entry[0]) for entry in n_most_common_strings]  # basis states
+    y_axis_values = [entry[1] for entry in n_most_common_strings]  # counts
+    num_successful_shots = sum(y_axis_values[:2])
+    print(f"Success ratio {num_successful_shots/n_shots} ")
+
+    fig = plt.figure()
+    ax = fig.add_axes([0, 0, 1.5, 1])
+    color_list = ["green"] * 2 + (["orange"] * (len(x_axis_values) - 2))
+    ax.bar(
+        x=x_axis_values,
+        height=y_axis_values,
+        color=color_list,
+    )
+    ax.set_title(label="Maxcut Results")
+    plt.ylim([0, 0.25 * n_shots])
+    plt.xlabel("Basis State")
+    plt.ylabel("Number of Shots")
+    plt.show()