Dirac-3S Max Clique Benchmarking

import numpy as np
import os
import time
from eqc_direct.client import EqcClient
from eqc_direct.utils import *
ip_address="172.18.15.110"
client = EqcClient(ip_address=ip_address) 

def read_dimacs_clq(file_path):
    """Read a DIMACS ascii .clq graph and return the adjacency matrix A.
    A[i, j] == 1 iff (i+1, j+1) is an edge (vertices in the file are 1-indexed).
    Lines: 'c ...' comments, 'p edge <n> <m>' header, 'e <u> <v>' edges.
    """
    n = None
    edges = []
    with open(file_path, 'r') as f:
        for line in f:
            parts = line.split()
            if not parts:
                continue
            tag = parts[0]
            if tag == 'c':
                continue
            elif tag == 'p':
                # p edge n m   (some files use 'p col n m' / 'p clq n m')
                n = int(parts[2])
                m = int(parts[3])
            elif tag == 'e':
                edges.append((int(parts[1]), int(parts[2])))
    if n is None:
        raise ValueError(f"no 'p' header found in {file_path}")
    A = np.zeros((n, n), dtype=np.uint8)
    for u, v in edges:
        if u == v:
            continue  # ignore self-loops
        A[u - 1, v - 1] = 1
        A[v - 1, u - 1] = 1  # undirected -> symmetric
    return A, n, m
def build_clique_hamiltonian(A):
    A = np.asarray(A, dtype=float)
    return -A

def clique_number_from_energy(energy, sum_constraint=1.0):
    R = float(sum_constraint)
    return 1.0 / (1.0 + energy / (R * R))
def extract_clique(x, A, tol=1e-6, scale=None):
    """Read a vertex set off the support of x and check it really is a clique.
    `tol` is RELATIVE: v is in the support when x[v] > tol * scale, with scale
    defaulting to max(x). 
    
    Returns (vertices, is_clique). Vertices are 0-indexed into A; add 1 to get
    the labels used in the DIMACS file.
    """
    x = np.asarray(x, dtype=float)
    if scale is None:
        scale = x.max() if x.size else 0.0
    if scale <= 0.0:
        return np.array([], dtype=int), False
    vertices = np.flatnonzero(x > tol * scale)
    k = len(vertices)
    if k == 0:
        return vertices, False
    sub = np.asarray(A, dtype=bool)[np.ix_(vertices, vertices)]
    # A clique on k vertices has k*(k-1) ones in its (symmetric, zero-diagonal) block.
    is_clique = bool(sub.sum() == k * (k - 1))
    return vertices, is_clique
def _extend_to_maximal(chosen, cand, Ab, x):
    """Grow `chosen` while candidates remain, highest weight first.
    """
    while True:
        idx = np.flatnonzero(cand)
        if idx.size == 0:
            return chosen
        w = x[idx]
        best = w.max()
        ties = idx[w >= best - 1e-12 * max(1.0, abs(best))]
        if ties.size > 1:
            deg = Ab[np.ix_(ties, idx)].sum(1)
            v = int(ties[np.argmax(deg)])
        else:
            v = int(ties[0])
        chosen.append(v)
        cand &= Ab[v]
        cand[v] = False
def greedy_clique_from_weights(x, A):
    """Repair step: greedily grow a genuine clique, taking vertices in order of x.
    """
    Ab = np.asarray(A, dtype=bool)
    x = np.asarray(x, dtype=float)
    cand = np.ones(Ab.shape[0], dtype=bool)
    return np.array(sorted(_extend_to_maximal([], cand, Ab, x)), dtype=int)

def direct_dirac3(indices, coefficients, relaxation_schedule, sum_constraint, num_samples, ip_address):
    eqc_client = EqcClient(ip_address=ip_address) #S9 "172.18.41.45", #S3 "172.18.41.173"
    lock_id, start_ts, end_ts=eqc_client.wait_for_lock()
    try:
        result_dict = eqc_client.solve_sum_constrained(
            lock_id=lock_id,
            poly_indices = indices,
            poly_coefficients = coefficients,
            relaxation_schedule = relaxation_schedule,
            sum_constraint = sum_constraint,
            num_samples = num_samples,
        )
        total_time = np.array(result_dict["preprocessing_time"])+np.array(result_dict["postprocessing_time"])+np.array(result_dict["runtime"])
        print(f"Total execution time(s):{total_time}")
    finally:
        # release lock when finished using the device
        lock_release_out = eqc_client.release_lock(lock_id=lock_id)
    return result_dict

instance_dir = "Instances/"
instance_name = "keller4.clq"
instance_path = os.path.join(instance_dir, instance_name)
try:
    A, n, m = read_dimacs_clq(instance_path)
    print("file loaded sucessfully")
except FileNotFoundError:
    print(f"{instance_path} does not exist")

file loaded sucessfully

name = os.path.splitext(instance_name)[0]
print(f"{name}: n = {n} vertices, m = {m} edges (declared)")
print(f"A shape: {A.shape}")
print(f"edges in A: {int(A.sum()) // 2}")
print(f"symmetric: {np.array_equal(A, A.T)}")
print(f"self-loops: {int(np.trace(A))}")
print(f"density: {A.sum() / (n * (n - 1)):.5f}")

keller4: n = 171 vertices, m = 9435 edges (declared)
A shape: (171, 171)
edges in A: 9435
symmetric: True
self-loops: 0
density: 0.64912

c = np.zeros(n)
H = build_clique_hamiltonian(A)
sum_constraint = 1
num_samples = 100
relaxation_schedule = 1

poly_indices, poly_coefficients = convert_hamiltonian_to_poly_format(
                    linear_terms=c,
                    quadratic_terms=H,   # the util mutates its input
                )
dirac_start = time.time()
print("Submitting job to Dirac-3S")
response = direct_dirac3(indices=poly_indices,
                        coefficients=poly_coefficients,
                        relaxation_schedule=relaxation_schedule,
                        sum_constraint=sum_constraint,
                        ip_address=ip_address,
                        num_samples=num_samples)
print(f"Dirac-3 Run complete. Response received.")
energies = response['energy']
solutions = response['solution']
runtimes = response['runtime']
preprocessing_time = response['preprocessing_time']
postprocessing_time = response['postprocessing_time']
print(f"Dirac-3 all samples energies:{[round(e, 2) for e in energies]}")
print(f"Dirac-3 all samples runtimes:{[round(t, 2) for t in runtimes]}")

WARNING:root:Max precision for EQC device is float32 input type was dtype float64. Input matrix will be rounded

Submitting job to Dirac-3S

/home/sutapa/Documents/eqc-direc-2.0.3/.venv/lib/python3.10/site-packages/eqc_direct/client.py:565: Warning: Max precision for EQC device is float32 input type was dtype float64. Input matrix will be rounded
  warnings.warn(warn_dtype_msg, Warning)

Total execution time(s):[2.62123575 2.58713802 2.01370343 2.7587886  2.04803766 2.67565884
 2.23382308 2.79223963 2.67314126 2.64704813 2.84181871 2.49983934
 2.79587569 2.90989884 2.82071676 2.66952636 2.52615025 1.86245462
 2.70607432 2.54609944 2.66275682 2.76443251 2.69481905 2.16523417
 2.80123534 2.7616445  2.03196877 2.04042424 2.80532928 2.7974691
 2.62519749 2.03304568 2.78892185 2.24321626 2.80168899 2.83994545
 2.62951382 2.03748437 2.54187473 2.63020476 2.86996313 2.68325921
 2.79529806 2.86349048 1.97953007 2.16283818 2.18016803 2.97041192
 2.82335903 2.07952028 2.68332769 2.79461773 2.54081376 2.40717832
 2.45432718 2.01539499 2.68938633 2.20637588 2.57783373 2.75900568
 2.50343429 2.69706388 2.68904926 2.06637112 2.11271995 2.66742732
 2.26051008 2.03095476 2.69906934 2.74248256 2.79104366 1.99890017
 2.78397586 2.7336476  2.65583034 2.90891715 2.6113561  2.27965211
 2.732641   2.80390247 2.15183202 2.81881728 2.88916387 2.75707759
 2.90165343 2.16630219 2.73159422 2.72164734 2.79233469 2.78964468
 2.03095688 1.87199597 2.63561304 2.66995545 2.76049509 2.82846047
 2.76028766 2.23443906 2.57353481 2.42797192]
Dirac-3 Run complete. Response received.
Dirac-3 all samples energies:[-0.86, -0.88, -0.88, -0.86, -0.86, -0.87, -0.86, -0.86, -0.86, -0.86, -0.86, -0.86, -0.86, -0.86, -0.86, -0.86, -0.87, -0.86, -0.87, -0.86, -0.86, -0.86, -0.86, -0.86, -0.86, -0.86, -0.86, -0.86, -0.86, -0.86, -0.86, -0.89, -0.86, -0.88, -0.86, -0.86, -0.88, -0.88, -0.86, -0.86, -0.86, -0.86, -0.86, -0.86, -0.86, -0.88, -0.86, -0.86, -0.86, -0.88, -0.86, -0.86, -0.88, -0.87, -0.87, -0.86, -0.86, -0.86, -0.86, -0.86, -0.87, -0.86, -0.86, -0.86, -0.86, -0.86, -0.86, -0.86, -0.86, -0.86, -0.86, -0.89, -0.86, -0.86, -0.86, -0.86, -0.87, -0.86, -0.86, -0.87, -0.88, -0.86, -0.86, -0.86, -0.86, -0.86, -0.86, -0.86, -0.86, -0.86, -0.88, -0.89, -0.86, -0.88, -0.86, -0.86, -0.86, -0.86, -0.86, -0.86]
Dirac-3 all samples runtimes:[1.97, 1.95, 2.0, 2.12, 2.04, 2.05, 2.06, 2.13, 2.03, 2.0, 2.19, 2.25, 2.16, 2.26, 2.18, 2.02, 1.89, 1.85, 2.07, 1.91, 2.01, 2.12, 2.05, 2.15, 2.15, 2.1, 2.02, 2.03, 2.15, 2.14, 1.98, 2.02, 2.14, 2.23, 2.14, 2.19, 1.99, 2.03, 1.89, 1.98, 2.22, 2.04, 2.16, 2.21, 1.93, 2.15, 2.17, 2.32, 2.19, 2.07, 2.05, 2.14, 1.91, 2.3, 1.82, 2.0, 2.04, 2.1, 1.92, 2.1, 2.16, 2.05, 2.03, 2.05, 2.1, 2.02, 2.0, 2.02, 2.06, 2.11, 2.15, 1.99, 2.13, 2.2, 2.01, 2.26, 1.97, 2.08, 2.22, 2.18, 2.14, 2.17, 2.24, 2.1, 2.26, 2.15, 2.07, 2.08, 2.13, 2.14, 2.02, 1.86, 1.98, 2.04, 2.12, 2.18, 2.11, 2.22, 1.93, 1.94]

best = min(zip(energies, solutions, runtimes), key=lambda x: x[0])
dirac_energy, best_solution, best_runtime = best
omega = clique_number_from_energy(dirac_energy, sum_constraint)
support, isclique = extract_clique(best_solution, A)
# Round every restart, not just the lowest-energy one.
cliques = [greedy_clique_from_weights(sol, A) for sol in solutions]
greedy = max(cliques, key=len)

print(f"number of samples:{num_samples}, relaxatrion schedule:{relaxation_schedule}")
print(f"best energy over {num_samples} samples:{dirac_energy:.6f}")
print(f"support size:{len(support)}, support is a clique:{isclique}")
print(f"clique vertices (1-indexed):{(greedy + 1).tolist()}")
print(f"{'energy':>16}{'omega_est':>12}{'clique':>9}{'time (s)':>12}")
print(f"{dirac_energy:>16.6f}{omega:>12.3f}{len(greedy):>9}{best_runtime:>12.3f}")
# Never report a set that is not actually a clique.
assert extract_clique(np.isin(np.arange(n), greedy).astype(float), A)[1], \
    "greedy repair returned a non-clique"

number of samples:100, relaxatrion schedule:1
best energy over 100 samples:-0.888889
support size:15, support is a clique:False
clique vertices (1-indexed):[6, 22, 30, 39, 49, 74, 79, 105, 112]
          energy   omega_est   clique    time (s)
       -0.888889       9.000        9       1.986