Source code for eqc_models.graph.maxcut

# (C) Quantum Computing Inc., 2024.
import networkx as nx
import numpy as np
from .base import TwoPartitionModel
[docs]
class MaxCutModel(TwoPartitionModel):
[docs]
    def decode(self, solution: np.ndarray) -> np.ndarray:
        """ Override the default decoding to use a the max cut metric to determine a solution """
        Gprime, solution = determine_solution(self.G, solution)
        cut_size = len(self.G.edges) - len(Gprime.edges)
        return solution
    @property
    def J(self) -> np.ndarray:
        return self.quad_objective
    @property
    def C(self) -> np.ndarray:
        return self.linear_objective
    
    @property
    def H(self):
        return self.linear_objective, self.quad_objective
[docs]
    def partition(self, solution):
        """ Return a dictionary with the partition number of each node """
        
        partition_num = {}
        for i, u in enumerate(self.node_map):
            if solution[i] == 0:
                partition_num[u] = 1
            else:
                partition_num[u] = 2
        return partition_num
[docs]
    def getCutSize(self, partition):
        cut_size = 0
        for u, v in self.G.edges:
            if partition[u]!=partition[v]:
                cut_size += 1
        return cut_size
[docs]
    def costFunction(self):
        """ 
        Parameters
        -------------
        
        None
        
        Returns
        --------------
        
        :C: linear operator (vector array of coefficients) for cost function
        :J: quadratic operator (N by N matrix array of coefficients ) for cost function
        
        """
        G = self.G
        self.node_map = list(G.nodes)
        variables = self.variables
        n = len(variables)
        self.upper_bound = np.ones((n,))
        
        J = np.zeros((n, n), dtype=np.float32)
        h = np.zeros((n, 1), dtype=np.float32)
        for u, v in G.edges:
            J[u, v] += 1
            J[v, u] += 1
            h[u, 0] -= 1
            h[v, 0] -= 1
        return h, J
def get_graph(n, d):
    """ Produce a repeatable graph with parameters n and d """
    seed = n * d
    return nx.random_graphs.random_regular_graph(d, n, seed)
def get_partition_graph(G, solution):
    """
    Build the partitioned graph, counting cut size 
    :parameters: G : nx.DiGraph, solution : np.ndarray
    :returns: nx.DiGraph, int
    
    """
    cut_size = 0
    Gprime = nx.DiGraph()
    Gprime.add_nodes_from(G.nodes)
    for i, j in G.edges:
        if solution[i] != solution[j]:
            cut_size += 1
        else:
            Gprime.add_edge(i, j)
    return Gprime, cut_size
def determine_solution(G, solution):
    """
    Use a simple bisection method to determine the binary solution. Uses
    the cut size as the metric.
    Returns the partitioned graph and solution.
    :parameters: G : nx.DiGraph, solution : np.ndarray
    :returns: nx.DiGraph, np.ndarray
    """
    solution = np.array(solution)
    test_vals = np.copy(solution)
    test_vals.sort()
    lower = 0
    upper = solution.shape[0] - 1
    best_cut_size = 0
    best_graph = G
    best_solution = None
    while upper > lower:
        middle = (upper + lower) // 2
        threshold = test_vals[middle]
        test_solution = (solution>=threshold).astype(np.int32)
        Gprime, cut_size = get_partition_graph(G, test_solution)
        if cut_size > best_cut_size:
            best_cut_size = cut_size
            lower = middle
            best_solution = test_solution
            best_graph = Gprime
        else:
            upper = middle
    return best_graph, best_solution
def get_maxcut_H(G, t):
    """ 
    Return a Hamiltonian representing the Maximum Cut Problem. Scale the problem using `t`.
    Automatically adds a slack qudit.
    
    """
    n = len(G.nodes)
    J = np.zeros(shape=(n+1, n+1), dtype=np.float32)
    h = np.zeros(shape=(n+1,1), dtype=np.float32)
    for u, v in G.edges:
        J[u, v] += 1
        J[v, u] += 1
        J[u, u] = 1
        J[v, v] = 1
        h[u] -= 1
        h[v] -= 1
    J *= 1/t**2
    h *= 1/t
    H = np.hstack([h, J])
    return H