Dirac-3S Planted Benchmarking

import numpy as np
import os
from eqc_direct.client import EqcClient
from eqc_direct.utils import *
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 loadQ(path, n):
    tri = np.load(path, mmap_mode="r")
    Q = np.empty((n,n), dtype=tri.dtype)
    iu = np.triu_indices(n)
    Q[iu] = tri
    Q.T[iu] = tri
    return Q

# Instance details
optimal_energy = 12000
sum_constraint = 100
num_var = 2000
k = 44
ub =10
seed =100
name = f"STQP_n_{num_var}_k_{k}_R_{sum_constraint}_seed_{seed}_ub_{ub}"
instance_path = os.path.join(f"Instances/{name}.npy")
c = np.zeros(num_var)
try:
    Q= loadQ(instance_path,num_var)
    print("loaded file sucessfully.")
except FileNotFoundError:
    print(f"File {instance_path} does not exist.")

loaded file sucessfully.

# solver parameters
relaxation_schedule = 1
num_samples = 10

# solve using Dirac-3S
poly_indices, poly_coefficients = convert_hamiltonian_to_poly_format(
                linear_terms=c,
                quadratic_terms=Q,
            )
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']

Submitting job to Dirac-3S

WARNING:root:Max precision for EQC device is float32 input type was dtype float64. Input matrix will be rounded
/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):[29.90855955 28.52003086 26.55735099 18.47528472 23.55796547 25.75906217
 18.2335756  15.7012367  18.30005711 26.98886002]
Dirac-3 Run complete. Response received.

print(f"Dirac-3S all samples energies:{[round(e, 2) for e in energies]}")
print(f"Dirac-3 all samples runtimes:{[round(t, 2) for t in runtimes]}")

Dirac-3S all samples energies:[12000.0, 13021.03, 12969.42, 13030.29, 13038.24, 13035.86, 13048.64, 13118.84, 13120.48, 13079.18]
Dirac-3 all samples runtimes:[27.94, 20.14, 18.12, 16.52, 21.61, 17.34, 16.28, 13.75, 16.34, 25.03]

best = min(zip(energies, solutions, runtimes), key=lambda x: x[0])
best_energy, best_solution, best_runtime = best
arr = np.array(best_solution)
indices = np.where(arr> 1e-6)[0]

print(f"samples:{num_samples}, Relaxation schedule:{relaxation_schedule}")
print(f"best energy over {num_samples} samples:{best_energy:.6f}")
print(f"support size:{len(indices)}")
print(f"time taken by best smaple {best_runtime}")

samples:10, Relaxation schedule:1
best energy over 10 samples:11999.999023
support size:44
time taken by best smaple 27.940969467163086

# compute relative gap
tol = 1e-7
abs_gap = best_energy-optimal_energy
relative_gap = (round(abs_gap,6)*100)/round(optimal_energy,6)

# print result
if abs_gap>tol:
    print(f"Dirac-3S solution is not optimal")
    print(f"Absolute Gap:{abs_gap}")
    print(f"Relative Gap(%):{abs(relative_gap)}")
else:
    print(f"Optimal solution found:{best_energy}.")

Optimal solution found:11999.9990234.