Dirac-3S QPLIB Benchmarking

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

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

def read_qplib(path):
    with open(path) as fh:
        toks = [t for t in (line.split("#")[0].strip() for line in fh) if t]
    pos = 0
    def nxt():
        nonlocal pos
        val = toks[pos]
        pos += 1
        return val
    def defaulted(length):
        """Read a `default value` / `number of non-defaults` / entries block."""
        arr = np.full(length, float(nxt()))
        for _ in range(int(nxt())):
            i, v = nxt().split()
            arr[int(i) - 1] = float(v)
        return arr
    name = nxt()
    probtype = nxt()
    sense = nxt().lower()
    n = int(nxt())
    m = int(nxt())
    # quadratic terms of the objective (each unordered pair listed once)
    Q = np.zeros((n, n), dtype=np.float32)
    for _ in range(int(nxt())):
        i, j, v = nxt().split()
        i, j, v = int(i) - 1, int(j) - 1, float(v)
        if i == j:
            Q[i, i] = v
        else:
            Q[i, j] = Q[j, i] = 0.5 * v   # split the pair coefficient over both triangles
    b = defaulted(n)            # linear terms of the objective
    obj_const = float(nxt())
    # linear terms of the constraints
    A = np.zeros((m, n))
    for _ in range(int(nxt())):
        k, i, v = nxt().split()
        A[int(k) - 1, int(i) - 1] = float(v)
    inf = float(nxt())
    lhs, rhs = defaulted(m), defaulted(m)
    lb, ub = defaulted(n), defaulted(n)
    return dict(name=name, type=probtype, sense=sense, n=n, m=m,
                Q=Q, b=b, obj_const=obj_const,
                A=A, lhs=lhs, rhs=rhs, lb=lb, ub=ub, inf=inf)

inst_dir = "Instances/"
inst_name = "QPLIB_2761.qplib"
inst = read_qplib(os.path.join(inst_dir, inst_name))
n, Q, b = inst["n"], inst["Q"], inst["b"]
M = 0.5 * Q
c = b.copy()
poly_indices, poly_coefficients = convert_hamiltonian_to_poly_format(
    linear_terms=c,
    quadratic_terms=M.copy(),   # the util mutates its input
)

def objective(x):
    """QPLIB objective, same convention as the poly submitted to Dirac-3."""
    x = np.asarray(x, dtype=float)
    return float(x @ M @ x + c @ x) + inst["obj_const"]
sum_constraint = float(inst["rhs"][0])   # sum_i x_i = 1
relaxation_schedule = 1
num_samples = 100
poly_indices, poly_coefficients = convert_hamiltonian_to_poly_format(
                    linear_terms=c,
                    quadratic_terms=M.copy(),   # 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]}")

Submitting job to Dirac-3S
Total execution time(s):[2.60927448 3.92017277 3.66367798 5.65877393 3.50433762 3.31685099
 4.9969642  2.65017603 2.94848334 3.96039353 4.76025672 5.69540859
 3.56813195 4.78678866 3.77307181 3.79501493 3.80302066 3.06004788
 3.61366897 3.79009187 3.36575746 2.88107067 3.71976098 4.5675675
 3.07030962 4.16581185 4.04716632 4.73760869 3.84494819 2.83097038
 3.98634813 4.67104459 3.73876965 3.85616669 2.68340677 4.18464988
 3.98370418 4.75931722 3.90973101 4.86702825 4.46201827 2.94292931
 4.27756047 3.44825469 3.7916764  2.74540909 3.15322383 4.18542743
 2.83082731 3.97312964 2.9065745  4.4457614  3.55981524 3.19108973
 3.93139399 3.77664767 3.94058345 2.63672721 3.83822486 3.5642059
 3.71930424 3.01746587 2.71147231 5.44961055 3.03422832 3.50443695
 3.14953573 4.63366989 3.4780005  4.37385515 4.02699331 4.06825119
 3.804987   2.96381695 3.53761173 3.83295924 3.04909765 3.63119259
 4.00034818 3.93047291 4.8892151  2.73325073 2.97676449 4.75357767
 3.2108481  3.2900839  3.74374173 3.93122042 3.67027874 5.30804687
 4.6108509  4.86014954 4.89530536 3.84987223 3.32743471 3.97914636
 3.97654084 4.17954905 3.96836448 5.34683647]
Dirac-3 Run complete. Response received.
Dirac-3 all samples energies:[0.09, 0.04, 0.06, 0.06, 0.09, 0.04, 0.04, 0.11, 0.09, 0.08, 0.07, 0.04, 0.13, 0.05, 0.11, 0.07, 0.07, 0.06, 0.0, 0.08, 0.05, 0.05, 0.06, 0.07, 0.06, 0.0, 0.0, 0.05, 0.11, 0.15, 0.07, 0.07, 0.11, 0.11, 0.14, 0.11, 0.04, 0.07, 0.05, 0.05, 0.0, 0.1, 0.0, 0.04, 0.05, 0.04, 0.09, 0.0, 0.1, 0.0, 0.07, 0.07, 0.07, 0.04, 0.11, 0.0, 0.05, 0.1, 0.04, 0.0, 0.05, 0.06, 0.07, 0.05, 0.09, 0.0, 0.05, 0.07, 0.0, 0.05, 0.06, 0.04, 0.05, 0.06, 0.02, 0.04, 0.05, 0.0, 0.01, 0.0, 0.04, 0.06, 0.06, 0.0, 0.07, 0.04, 0.04, 0.09, 0.04, 0.0, 0.0, 0.0, 0.04, 0.0, 0.06, 0.05, 0.0, 0.07, 0.0, 0.05]
Dirac-3 all samples runtimes:[2.38, 3.69, 3.43, 5.43, 3.27, 3.09, 4.77, 2.42, 2.72, 3.73, 4.53, 5.47, 3.34, 4.56, 3.54, 3.57, 3.57, 2.83, 3.38, 3.56, 3.14, 2.65, 3.49, 4.34, 2.84, 3.94, 3.82, 4.51, 3.61, 2.6, 3.76, 4.44, 3.51, 3.63, 2.45, 3.95, 3.75, 4.53, 3.68, 4.64, 4.23, 2.71, 4.05, 3.22, 3.56, 2.52, 2.92, 3.96, 2.6, 3.74, 2.68, 2.63, 3.33, 2.96, 3.7, 3.55, 3.71, 2.41, 3.61, 3.33, 3.49, 2.79, 2.48, 3.64, 2.8, 3.27, 2.92, 4.4, 3.25, 4.14, 3.8, 3.84, 3.58, 2.73, 3.31, 3.6, 2.82, 3.4, 3.77, 3.7, 4.66, 2.5, 2.75, 4.52, 2.98, 3.06, 3.51, 3.7, 3.44, 5.08, 4.38, 4.63, 4.67, 3.62, 3.1, 3.75, 3.75, 3.95, 3.74, 5.12]

best = min(zip(energies, solutions, runtimes), key=lambda x: x[0])
dirac_energy, best_solution, best_runtime = best
print(f"Eenergy value obtained by Dirac-3:{dirac_energy}")
print(f"Solution sum:{sum(best_solution)} (target {sum_constraint}), support:{np.count_nonzero(best_solution)}")
print(f"Preprocessing Time:{preprocessing_time}")
print(f"Time taken for this run:{best_runtime}")

Eenergy value obtained by Dirac-3:0.0010485
Solution sum:1.0000001 (target 1.0), support:2
Preprocessing Time:0.22956259548664093
Time taken for this run:3.384004592895508