Solve Max-Cut using Fire Opal

import fireopal as fo
import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
from pyscipopt import Model
from sympy import symbols, Poly
import qctrlvisualizer as qv

# These are the properties for the publicly available provider for IBM backends.
# If you have access to a private provider and wish to use it, replace these values.
hub = "ibm-q"
group = "open"
project = "main"
token = "YOUR_IBM_TOKEN"
credentials = fo.credentials.make_credentials_for_ibmq(
    token=token, hub=hub, group=group, project=project
)

# Enter your desired IBM backend here.
backend_name = "desired_backend"

def maxcut_obj(bitstring: str, graph: nx.Graph) -> float:
    """
    Given a bitstring, this function returns the number of edges
    shared between the two partitions of the graph.
    """
    obj = 0
    # Iterate through each edge in the graph
    for i, j, data in graph.edges(data=True):
        # Calculate the contribution of the edge to the Max-Cut objective
        # If the nodes are in different partitions, the term (int(bitstring[i]) - int(bitstring[j]))**2 will be 1
        # Otherwise, it will be 0
        obj -= data.get("weight", 1) * (int(bitstring[i]) - int(bitstring[j])) ** 2
    return obj


def get_cost_distribution(
    distribution: dict[str, int], graph: nx.Graph
) -> tuple[np.ndarray, np.ndarray]:
    """
    Given a bitstring count distribution and a graph, returns the
    costs and relative probabilities as numpy arrays.
    """
    # Total number of bitstrings (shots)
    shot_count = sum(distribution.values())
    costs = []
    probabilities = []
    # Calculate the cost and probability for each bitstring
    for bitstring, count in distribution.items():
        # Reverse the bitstring to match graph node indices
        cost = maxcut_obj(bitstring[::-1], graph)
        costs.append(cost)
        probabilities.append(count / shot_count)
    return np.array(costs), np.array(probabilities)


def generate_random_bitstrings(length, num_bitstrings, rng_seed=0) -> list[str]:
    """
    Generate a random sampling of bitstrings.
    """
    # Initialize the random number generator with a seed for reproducibility
    rng = np.random.default_rng(seed=rng_seed)
    # Generate random bitstrings
    random_array = rng.integers(0, 2, size=(num_bitstrings, length))
    bitstrings = ["".join(row.astype(str)) for row in random_array]
    return bitstrings


def calculate_percent_gap_to_optimal(
    costs: np.ndarray, optimal_cut_value: float
) -> None:
    """
    Determine the percent gap between the optimal cut value and the best cut value
    found from Fire Opal.
    """
    # Find the best cut value found (most negative cost)
    best_cut_value_found = max(-costs)
    # Calculate the percentage difference from the optimal value
    percent_optimality = np.round(
        100 - (abs(optimal_cut_value - best_cut_value_found) / optimal_cut_value * 100),
        2,
    )
    # Print the results
    print(
        f"The best cut value found by Fire Opal is {best_cut_value_found}. \n"
        f"This solution is {percent_optimality}% optimized relative to the classically determined solution, {optimal_cut_value}."
    )


def analyze_maxcut_results(
    distribution: dict[str, int], graph: nx.Graph, optimal_cut_value: float
) -> None:
    """
    Analyze results from Max-Cut execution.
    """
    # Number of nodes in the graph
    qubit_count = graph.number_of_nodes()
    # Total number of bitstrings (shots)
    shot_count = sum(distribution.values())

    # Get cost distribution from the results
    costs, probs = get_cost_distribution(distribution, graph)
    # Generate random bitstrings and get their cost distribution
    random_bitstrings = generate_random_bitstrings(qubit_count, shot_count)
    random_probs = np.ones(shot_count) / shot_count
    random_costs = np.asarray(
        [maxcut_obj(bitstring[::-1], graph) for bitstring in random_bitstrings]
    )

    # Plot histograms of cut values
    plt.style.use(qv.get_qctrl_style())
    fig, ax = plt.subplots(1, 1)

    n, bin_edges, _ = ax.hist(
        [-costs, -random_costs],
        np.arange(-max(random_costs), -min(costs) + 1, 1),
        weights=[probs, random_probs],
        label=["Fire Opal QAOA solver", "Random Sampling"],
        color=[qv.QCTRL_STYLE_COLORS[0], qv.QCTRL_STYLE_COLORS[1]],
    )
    ax.set_xlabel("Cut Value")
    ax.set_ylabel("Probability")
    ax.legend()

    # Calculate and print the percent gap to optimal
    calculate_percent_gap_to_optimal(costs=costs, optimal_cut_value=optimal_cut_value)


def solve_maxcut_classically(graph):
    """
    Solve the Max-Cut problem classically using linear programming.
    """
    model = Model("MaxCut")
    # Add binary variables for edges
    e = {(u, v): model.addVar(vtype="B", name=f"e({u, v})") for u, v in graph.edges()}
    # Add binary variables for nodes
    x = {u: model.addVar(vtype="B", name=f"x({u})") for u in graph.nodes()}

    # Define the objective function to maximize the cut value
    objective = sum(
        data.get("weight", 1) * e[(u, v)] for u, v, data in graph.edges(data=True)
    )
    model.setObjective(objective, "maximize")

    # Add constraints to ensure valid cuts
    for u, v in graph.edges():
        model.addCons(e[(u, v)] <= x[u] + x[v])
        model.addCons(e[(u, v)] <= 2 - (x[u] + x[v]))

    model.hideOutput(True)
    model.optimize()

    # Ensure the solution is optimal
    assert model.getStatus() == "optimal"

    # Return the solution indicating which partition each node belongs to
    return {u: model.getVal(x[u]) for u in graph.nodes()}


def cut_value(graph, cut):
    """
    Calculate the cut value of a given partition.
    """
    value = 0
    # Calculate the total weight of the cut edges
    for i, j, data in graph.edges(data=True):
        value += data.get("weight", 1) * np.round(cut[i] - cut[j]) ** 2
    return value

import networkx as nx
import numpy as np

# Generate a random graph with 156 nodes
unweighted_graph = nx.random_regular_graph(d=3, n=156, seed=8)

# Optionally, visualize the graph
nx.draw(unweighted_graph, nx.kamada_kawai_layout(unweighted_graph), node_size=100)

unweighted_job = fo.solve_qaoa(
    problem=unweighted_graph,
    credentials=credentials,
    problem_type="maxcut",
    backend_name=backend_name,
)

unweighted_result = unweighted_job.result()

solution = solve_maxcut_classically(unweighted_graph)
optimal_cut_value = cut_value(unweighted_graph, solution)
print(f"Classically derived cut value: {optimal_cut_value}")

analyze_maxcut_results(
    distribution=unweighted_result["final_bitstring_distribution"],
    graph=unweighted_graph,
    optimal_cut_value=optimal_cut_value,
)

from fireopal import print_package_versions

print_package_versions()