graph utils and bug fixes (#27)

* minimum edge cover

* graph utils

* Grover for 1 qubit

* do not use end point

* bugfixes

Author: fgfuchs

examples/plotroutines.py

@@ -1,4 +1,4 @@
-import matplotlib.pyplot as pl
+import matplotlib.pyplot as plt
 from mpl_toolkits.axes_grid1 import make_axes_locatable
 from matplotlib.ticker import MaxNLocator
 
@@ -11,15 +11,15 @@
 
 def __plot_landscape(A, extent, fig):
     if not fig:
-        fig = pl.figure(figsize=(6, 6), dpi=80, facecolor="w", edgecolor="k")
-    _ = pl.xlabel(r"$\gamma$")
-    _ = pl.ylabel(r"$\beta$")
+        fig = plt.figure(figsize=(6, 6), dpi=80, facecolor="w", edgecolor="k")
+    _ = plt.xlabel(r"$\gamma$")
+    _ = plt.ylabel(r"$\beta$")
     ax = fig.gca()
-    _ = pl.title("Expectation value")
+    _ = plt.title("Expectation value")
     im = ax.imshow(A, interpolation="bicubic", origin="lower", extent=extent)
     divider = make_axes_locatable(ax)
     cax = divider.append_axes("right", size="5%", pad=0.05)
-    _ = pl.colorbar(im, cax=cax)
+    _ = plt.colorbar(im, cax=cax)
 
 
 def plot_E(qaoa_instance, fig=None):
@@ -57,21 +57,21 @@ def plot_ApproximationRatio(
         exp = np.array(exp)
 
     if not fig:
-        ax = pl.figure().gca()
+        ax = plt.figure().gca()
     else:
         ax = fig.gca()
-    pl.hlines(1, 1, maxdepth, linestyles="solid", colors="black")
-    pl.plot(
+    plt.hlines(1, 1, maxdepth, linestyles="solid", colors="black")
+    plt.plot(
         np.arange(1, maxdepth + 1),
         (maxcost - exp) / (maxcost - mincost),
         style,
         label=label,
     )
-    pl.ylim(0, 1.01)
-    pl.xlim(1 - 0.25, maxdepth + 0.25)
-    _ = pl.ylabel("appr. ratio")
-    _ = pl.xlabel("depth")
-    _ = pl.legend(loc="lower right", framealpha=1)
+    plt.ylim(0, 1.01)
+    plt.xlim(1 - 0.25, maxdepth + 0.25)
+    _ = plt.ylabel("appr. ratio")
+    _ = plt.xlabel("depth")
+    _ = plt.legend(loc="lower right", framealpha=1)
     ax.xaxis.set_major_locator(MaxNLocator(integer=True))
 
 
@@ -83,21 +83,21 @@ def plot_successprob(qaoa_instance, maxdepth, label, style="", fig=None, shots=1
     successp = np.array(successp)
 
     if not fig:
-        ax = pl.figure().gca()
+        ax = plt.figure().gca()
     else:
         ax = fig.gca()
-    pl.hlines(1, 1, maxdepth, linestyles="solid", colors="black")
-    pl.plot(
+    plt.hlines(1, 1, maxdepth, linestyles="solid", colors="black")
+    plt.plot(
         np.arange(1, maxdepth + 1),
         successp,
         style,
         label=label,
     )
-    pl.ylim(0, 1.01)
-    pl.xlim(1 - 0.25, maxdepth + 0.25)
-    _ = pl.ylabel("success prob")
-    _ = pl.xlabel("depth")
-    _ = pl.legend(loc="lower right", framealpha=1)
+    plt.ylim(0, 1.01)
+    plt.xlim(1 - 0.25, maxdepth + 0.25)
+    _ = plt.ylabel("success prob")
+    _ = plt.xlabel("depth")
+    _ = plt.legend(loc="lower right", framealpha=1)
     ax.xaxis.set_major_locator(MaxNLocator(integer=True))
 
 
@@ -126,24 +126,62 @@ def plot_angles(qaoa_instance, depth, label, style="", fig=None):
     angles = qaoa_instance.optimization_results[depth].get_best_angles()
 
     if not fig:
-        ax = pl.figure().gca()
+        ax = plt.figure().gca()
     else:
         ax = fig.gca()
 
-    pl.plot(
+    plt.plot(
         np.arange(1, depth + 1),
         angles[::2],
         "--" + style,
         label=r"$\gamma$ " + label,
     )
-    pl.plot(
+    plt.plot(
         np.arange(1, depth + 1),
         angles[1::2],
         "-" + style,
         label=r"$\beta$ " + label,
     )
-    pl.xlim(1 - 0.25, depth + 0.25)
-    _ = pl.ylabel("parameter")
-    _ = pl.xlabel("depth")
-    _ = pl.legend()
+    plt.xlim(1 - 0.25, depth + 0.25)
+    _ = plt.ylabel("parameter")
+    _ = plt.xlabel("depth")
+    _ = plt.legend()
     ax.xaxis.set_major_locator(MaxNLocator(integer=True))
+
+
+def draw_colored_graph(G, edge_colors):
+    # Draw the graph with colored edges
+    # extend the color_map if necessary
+    color_map = [
+        "red",
+        "blue",
+        "green",
+        "purple",
+        "orange",
+        "pink",
+        "brown",
+        "gray",
+        "yellow",
+        "cyan",
+    ]
+    pos = nx.spring_layout(G)  # Positions for all nodes
+
+    # Draw nodes
+    nx.draw_networkx_nodes(G, pos, node_size=700, node_color="lightgray")
+
+    # Draw edges with colors
+    for color_idx, edges in edge_colors.items():
+        nx.draw_networkx_edges(
+            G,
+            pos,
+            edgelist=edges,
+            width=2,
+            edge_color=color_map[color_idx % len(color_map)],
+        )
+
+    # Draw labels
+    nx.draw_networkx_labels(G, pos, font_size=20, font_family="sans-serif")
+
+    # Show the graph
+    plt.axis("off")
+    plt.show()

qaoa/mixers/grover_mixer.py

@@ -13,6 +13,7 @@ def __init__(self, subcircuit: InitialState) -> None:
         """
         self.subcircuit = subcircuit
         self.mixer_param = Parameter("x_beta")
+        self.N_qubits = subcircuit.N_qubits
 
     def create_circuit(self):
         # given feasibel states f \in F,
@@ -27,7 +28,12 @@ def create_circuit(self):
         # X^n
         self.circuit.x(range(self.subcircuit.N_qubits))
         # C^{n-1}Phase
-        phase_gate = PhaseGate(-self.mixer_param).control(self.subcircuit.N_qubits - 1)
+        if self.subcircuit.N_qubits == 1:
+            phase_gate = PhaseGate(-self.mixer_param)
+        else:
+            phase_gate = PhaseGate(-self.mixer_param).control(
+                self.subcircuit.N_qubits - 1
+            )
         self.circuit.append(phase_gate, self.circuit.qubits)
         # X^n
         self.circuit.x(range(self.subcircuit.N_qubits))

qaoa/problems/qubo_problem.py

@@ -28,6 +28,8 @@ def __init__(self, Q=None, c=None, b=None) -> None:
         )
         n = Q.shape[0]
 
+        self.N_qubits = n
+
         # Check if Q is lower triangular
         self.lower_triangular_Q = np.allclose(Q, np.tril(Q))
 

qaoa/qaoa.py

@@ -288,9 +288,9 @@ def sample_cost_landscape(
         depth = 1
 
         tmp = angles["gamma"]
-        self.gamma_grid = np.linspace(tmp[0], tmp[1], tmp[2])
+        self.gamma_grid = np.linspace(tmp[0], tmp[1], tmp[2], endpoint=False)
         tmp = angles["beta"]
-        self.beta_grid = np.linspace(tmp[0], tmp[1], tmp[2])
+        self.beta_grid = np.linspace(tmp[0], tmp[1], tmp[2], endpoint=False)
 
         if self.backend.configuration().local:
             if self.sequential:

qaoa/util/__init__.py

@@ -1,3 +1,4 @@
 from .statistic import Statistic
 from .flip import BitFlip
 from .post import *
+from .graphutils import GraphHandler

qaoa/util/flip.py

@@ -7,7 +7,7 @@ def __init__(self, n):
         self.circuit = None
         self.N_qubits = n
 
-    def boost_samples(self, problem, string, K = 5):
+    def boost_samples(self, problem, string, K=5):
         """
         Random bitflips on string/list of strings to increase cost.
 
@@ -22,13 +22,13 @@ def boost_samples(self, problem, string, K = 5):
         old_string = string
         cost = problem.cost(string[::-1])
 
-        for _ in range (K):
+        for _ in range(K):
             shuffled_indices = np.arange(self.N_qubits)
             np.random.shuffle(shuffled_indices)
 
             for i in shuffled_indices:
                 string_arr_altered = np.copy(string_arr)
-                string_arr_altered[i] = not(string_arr[i])
+                string_arr_altered[i] = not (string_arr[i])
                 string_altered = "".join(map(str, string_arr_altered))
                 new_cost = problem.cost(string_altered[::-1])
 

qaoa/util/graphutils.py

@@ -0,0 +1,110 @@
+import networkx as nx
+
+
+class GraphHandler:
+    def __init__(self, G):
+        if isinstance(G, nx.DiGraph):
+            raise Exception("Graph should be undirected.")
+        if any(degree <= 1 for node, degree in G.degree()):
+            print(
+                "Graph contains nodes with one or zero edges. These can be removed to reduce the size of the problem."
+            )
+
+        self.num_nodes = G.number_of_nodes()
+        self.num_edges = G.number_of_edges()
+
+        # ensure graph has labels 0, 1, ..., num_V-1
+        G_int = self.__ensure_integer_labels__(G)
+        # relabel to make node n-1 the one with maximum degree
+        self.G = self.__get_graph_maxdegree_last_node__(G_int)
+        # to avoid a deep circuit, we partition the edges into sets which can be executed in parallel
+        if not nx.is_isomorphic(G, self.G):
+            raise Exception("Something went wrong.")
+
+        self.__minimum_edge_coloring__()
+
+    def __ensure_integer_labels__(self, G):
+        # Check if nodes are already labeled as 0 to num_nodes-1
+        if set(G.nodes) == set(range(self.num_nodes)):
+            return (
+                G  # Return the graph unchanged if nodes are already labeled correctly
+            )
+
+        # If nodes are not labeled correctly, create a mapping to relabel them
+        node_mapping = {node: i for i, node in enumerate(G.nodes)}
+
+        # Create a new graph with relabeled nodes
+        H = nx.relabel_nodes(G, node_mapping, copy=True)
+
+        return H
+
+    def __map_colors_to_edges__(self, line_graph_colors, original_graph):
+        # Map colors to edges in the original graph G
+        color_to_edges = {}
+
+        # Each node in the line graph corresponds to an edge in the original graph G
+        for edge_in_line_graph, color in line_graph_colors.items():
+            original_edge = edge_in_line_graph  # This is the corresponding edge in G
+            if color not in color_to_edges:
+                color_to_edges[color] = []
+            color_to_edges[color].append(original_edge)
+
+        # Perform consistency check
+
+        # Get the set of all edges in G
+        edges_in_G = set(self.G.edges())
+        # Get the set of all edges in color_to_edges
+        edges_in_coloring = set(
+            edge for edges in color_to_edges.values() for edge in edges
+        )
+
+        if edges_in_G != edges_in_coloring:
+            raise ValueError(
+                "The colored edges do not match the edges in the original graph!"
+            )
+
+        return color_to_edges
+
+    def __get_graph_maxdegree_last_node__(self, G):
+        # Get node of highest degree
+        j = sorted(G.degree(), key=lambda x: x[1], reverse=True)[0][0]
+        if j == self.num_nodes - 1:
+            return G
+        else:
+            # Create a mapping to swap node j and n-1
+            mapping = {j: self.num_nodes - 1, self.num_nodes - 1: j}
+
+            # Relabel the nodes
+            H = nx.relabel_nodes(G, mapping, copy=True)
+
+            return H
+
+    def __minimum_edge_coloring__(self, repetitions=100):
+        #
+        # a graph G
+        # returns minimum edge coloring, i.e., a dict containting the edges for each color
+        # this can be used to minimize the depth needed to implement diagonal cost Hamiltonians
+        # example output
+        # { 3: [(0, 1), (2, 7), (3, 9), (5, 6)],
+        #   1: [(0, 3), (1, 5), (7, 9)],
+        #   2: [(0, 9), (1, 6), (2, 5), (3, 8), (4, 7)],
+        #   0: [(1, 4), (2, 9), (3, 7), (5, 8)] }
+        #
+
+        # Convert the graph to its line graph
+        line_G = nx.line_graph(self.G)
+
+        ncolors = self.num_edges + 1
+        for _ in range(repetitions):
+            # Apply greedy vertex coloring on the line graph
+            line_graph_colors = nx.coloring.greedy_color(
+                line_G, strategy="random_sequential"
+            )
+
+            # groups of parallel edges
+            pe = self.__map_colors_to_edges__(line_graph_colors, self.G)
+
+            num_classes = len(pe)
+            if num_classes < ncolors:
+                self.parallel_edges = pe
+                ncolors = num_classes

qaoa/util/post.py

@@ -1,21 +1,22 @@
 import statistics as stat
 import numpy as np
 
+
 def post_processing(instance, samples, K=5):
     """Performs classical post-processing on bitstrings by applying random bit flips.
     Resets and updates self.stat
 
     input:
         instance: instance to perform function on
-        samples (dict or list or str): The bitstring(s) to be processed. 
+        samples (dict or list or str): The bitstring(s) to be processed.
             The acceptable formats are:
             - dict{bitstring: count}: A dictionary with bitstrings as keys and their counts as values.
             - list(bitstring, bitstring, ...): A list of bitstrings.
             - str(bitstring): A single bitstring.
         K (int): The number of times to iterate through each bitstring and apply random bit flips.
 
     returns:
-        dict: A dictionary with the altered bitstrings as keys and their counts as values. 
+        dict: A dictionary with the altered bitstrings as keys and their counts as values.
         If no better bitstring is found, the original bitstring is the key.
     """
     instance.stat.reset()
@@ -26,19 +27,20 @@ def post_processing(instance, samples, K=5):
 
     for string in samples:
         boosted = instance.flipper.boost_samples(
-            problem=instance.problem,
-            string=string,
-            K=K            
+            problem=instance.problem, string=string, K=K
         )
         try:
             count = samples[string]
         except:
             count = 1
 
-        instance.stat.add_sample(instance.problem.cost(boosted[::-1]), count, boosted[::-1])
+        instance.stat.add_sample(
+            instance.problem.cost(boosted[::-1]), count, boosted[::-1]
+        )
         hist_post[boosted] = hist_post.get(boosted, 0) + count
     return hist_post
 
+
 def post_process_all_depths(instance, K=5):
     """Performs post-processing of job.result().get_counts() 100 times after each layer.
 
@@ -56,7 +58,7 @@ def post_process_all_depths(instance, K=5):
         for i in range(100):
             post_processing(
                 instance=instance,
-                samples=hist, 
+                samples=hist,
                 K=K,
             )
             exp_in_layers[d] = exp_in_layers.get(d, []) + [-instance.stat.get_CVaR()]