Team Pineapple Submission

team_solutions/pineapple/Final.py

@@ -0,0 +1,72 @@
+from matplotlib import pyplot as plt
+from pytket.extensions.qiskit import AerBackend
+from pytket.backends.backend import Backend
+from pytket.backends.backendresult import BackendResult
+from pytket.passes import DecomposeBoxes
+from pytket.utils import gen_term_sequence_circuit
+import numpy as np
+from pytket import Qubit, Circuit
+from pytket.pauli import QubitPauliString, Pauli
+from pytket.utils import QubitPauliOperator
+from typing import List, Tuple, Callable
+import networkx as nx
+import networkx.algorithms.isomorphism.vf2userfunc as vf2
+import networkx as nx
+import numpy as np
+import matplotlib.pyplot as plt
+import networkx.algorithms.isomorphism.vf2userfunc as vf2
+import os
+import pickle
+
+
+def subgraphInduce(G, edge, depth, rename=True):
+    # return: subgraph induced by edge
+    current = set(edge)
+    edges = set([edge])
+    for i in range(depth):
+        next = set()
+        for node in current:
+            next.update(G.neighbors(node))
+            edges.update(G.edges(node))
+        current.update(next)
+    if rename:
+        # Rename nodes to 0, 1, 2, ... such that 0 and 1 are the central edge.
+        current.remove(edge[0])
+        current.remove(edge[1])
+        current = [edge[0], edge[1]] + list(current)
+
+        edges = [(current.index(e[0]), current.index(e[1]))
+                 for e in edges]
+    return nx.Graph(list(edges))
+
+
+x = pickle.load(open("subgraphs3.pkl", "rb"))  # Subgraphs with frequency
+
+our_subgraphs = pickle.load(open("subgraphs.pkl", "rb"))  # Subgraphs
+
+print(len(our_subgraphs))
+freq = []
+for idx, subgraph in enumerate(our_subgraphs):
+    found = False
+    for i in x:
+        if nx.is_isomorphic(i[0], nx.Graph(subgraph)):
+            found = True
+            freq.append([idx, i[1]])
+            break
+    if not found:
+        freq.append([idx, 0])
+
+
+def filename(graph):
+    h = "".join(str(i[0]) + str(i[1]) for i in graph)
+    return h
+
+
+print(len(freq))
+freq = sorted(freq, key=lambda x: x[1], reverse=True)
+print("\n".join([str(i) for i in freq]))
+
+# Figure out which ones we urgently need to find
+for i, f in freq:
+    if not os.path.exists(str(filename(our_subgraphs[i])) + ".pkl"):
+        print("FUCK", i, f)

team_solutions/pineapple/README.md

@@ -0,0 +1,21 @@
+## Challenges on QAOA 
+
+1. This QAOA [(Quantum Approximate Optimisation algorithm)](https://arxiv.org/abs/1411.4028) implementaion uses the most naive possible classical optimisation strategy. Parameters are sampled from a uniform distribution and if a list of parameters increases the value of the cost function these values are stored as the best guess so far. Can you improve on this using a more sophisticated optimisation strategy? COBAYLA and SPSA are two possible methods.
+
+2. The maxcut problem is one very common application of QAOA. Can you create an implementation of QAOA applied to a different problem? Examples of such problems included 3SAT and the maximum clique problem. Perhaps try and create and implementation which works for a more general Hamiltonian that could contain non-commuting Pauli terms like those found in Quantum Chemistry. Think about what additonal complexity would be added by a Hamiltonian with non-commuting terms. Interesting Hamiltonians to consider could be the Transverse Field Ising Model (TFIM), diatomic Hydrogen or a simple compound like lithium hydride. 
+
+3. The given code implements QAOA on the idealised AerBackend simulator. Try instead to use a device/emulator with noise (i.e. the H1-2 emulator with the pytket-quantinuum extension). Can you optimise your circuit with pytket passes to improve performance in the presence of noise?
+
+4. Currently the circuits have to be recompiled on every iteration leading to a non-trivial compilation overhead if we use a large number of iterations. Can you think of a way to improve this? 
+
+5. Implement a qubit reuse strategy to allow for the execution of a large QAOA instance on a small quantum device/emulator. See [this paper](https://arxiv.org/abs/2210.08039) for ideas.
+
+## Resources
+
+1. QAOA original paper (Farhi et al) -> https://arxiv.org/abs/1411.4028
+2. Quantinuum Qubit reuse paper (DeCross et al) -> https://arxiv.org/abs/2210.08039
+3. pytket API documentation -> https://cqcl.github.io/tket/pytket/api/
+4. User manual -> https://cqcl.github.io/pytket/manual/index.html
+5. Notebook examples -> https://github.com/CQCL/pytket/tree/main/examples
+6. Qiskit textbook section on QAOA -> https://qiskit.org/textbook/ch-applications/qaoa.html
+7. Recent QAOA/qiskit review -> https://arxiv.org/abs/2301.09535

team_solutions/pineapple/Subgraph.py

@@ -0,0 +1,75 @@
+import networkx as nx
+import numpy as np
+import matplotlib.pyplot as plt
+import networkx.algorithms.isomorphism.vf2userfunc as vf2
+
+
+def subgraphInduce(G, edge, depth, rename=False):
+    # return: subgraph induced by edge
+    current = set([edge[0], edge[1]])
+    edges = set([edge])
+    for i in range(depth):
+        next = set()
+        for node in current:
+            next.update(G.neighbors(node))
+            edges.update(G.edges(node))
+        current.update(next)
+    if rename:
+        current = {x: i for i, x in enumerate(list(current))}
+        edges = [(current[edge[0]], current[edge[1]]) for edge in edges]
+    return nx.Graph(list(edges))
+
+
+# adj = list(G.neighbors(0))
+# print(adj[0])
+# print(subgraphInduce(G, (0, adj[0]), 1).edges())
+# print(subgraphInduce(G, (0, adj[0]), 1, True).edges())
+
+'''
+for i in range(40):
+    for t in range(20, 30, 2):
+        G = nx.random_regular_graph(3, t)
+        for edge in G.edges():
+            subgraph = subgraphInduce(G, edge, 2)
+            seen = False
+            # for k, (i, j) in enumerate(res):
+            #     if nx.is_isomorphic(subgraph, i):
+            #         seen = True
+            #         res[k][1] += 1
+            #         break
+            for i in res:
+                if nx.is_isomorphic(subgraph, i):
+                    seen = True
+                    break
+            if not seen:
+                #print(len(res))
+                #res.append([subgraph, 1])
+                res.append(subgraph)'''
+
+
+def gimme_subgraphs():
+    res = []
+
+    for t in range(20, 100, 2):
+        G = nx.random_regular_graph(3, t)
+        for edge in G.edges():
+            subgraph = subgraphInduce(G, edge, 2)
+            seen = False
+            for i in res:
+                if nx.is_isomorphic(subgraph, i):
+                    seen = True
+                    # res[i][1] += 1
+                    break
+            if not seen:
+                # print(len(res))
+                # res.append([subgraph, 1])
+                res.append(subgraph)
+    print("Generated Subgraphs len = ", len(res))
+    return res
+
+
+# res = sorted(res, key=lambda x: x[1], reverse=True)
+# for i in range(100):
+#    print(res[i][1])
+# nx.draw(G)
+# plt.show()

team_solutions/pineapple/interpolate_maxcut.py

@@ -0,0 +1,202 @@
+from maxcut_plotting import plot_maxcut_results
+from pytket.extensions.qiskit import AerBackend
+from pytket.backends.backend import Backend
+from pytket.backends.backendresult import BackendResult
+from pytket.passes import DecomposeBoxes
+from pytket.utils import gen_term_sequence_circuit
+import numpy as np
+from pytket import Qubit, Circuit
+from pytket.pauli import QubitPauliString, Pauli
+from pytket.utils import QubitPauliOperator
+from typing import List, Tuple, Callable
+import networkx as nx
+
+# Define graph.
+
+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)]
+
+
+def qaoa_graph_to_cost_hamiltonian(
+    edges: List[Tuple[int, int]], cost_angle: float
+) -> QubitPauliOperator:
+    """
+    This function takes a list of edges and a cost angle and returns a QubitPauliOperator
+    representing the cost Hamiltonian for the QAOA algorithm.
+
+    """
+    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)
+
+
+cost_angle = 1.0
+cost_ham_qpo = qaoa_graph_to_cost_hamiltonian(max_cut_graph_edges, cost_angle)
+print(cost_ham_qpo)
+
+
+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:
+    """
+    Create a QAOA circuit for the MaxCut problem.
+    """
+
+    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, maximize=False) -> float:
+    energy = 0.0
+    dist = results.get_distribution()
+    if maximize:
+        for meas in dist.keys():
+            energy = max(energy, sum((meas[i] ^ meas[j]) for i, j in edges))
+    else:
+        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
+
+
+backend = AerBackend()
+comp = backend.get_compiled_circuit
+
+res = qaoa_calculate(
+    backend,
+    backend.default_compilation_pass(2).apply,
+    shots=5000,
+    iterations=100,
+    seed=12345,
+    maximize=True
+)
+
+
+plot_maxcut_results(res, 6)