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
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

# 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()