{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "87b3ca42",
   "metadata": {},
   "source": [
    "# Dirac-3S Max Clique Benchmarking"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80e94f29",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import os\n",
    "import time\n",
    "from eqc_direct.client import EqcClient\n",
    "from eqc_direct.utils import *\n",
    "\n",
    "ip_address=\"172.18.15.110\"\n",
    "client = EqcClient(ip_address=ip_address) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26cb1b1d",
   "metadata": {},
   "source": [
    "### Read data file and build objective functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c16773c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_dimacs_clq(file_path):\n",
    "    \"\"\"Read a DIMACS ascii .clq graph and return the adjacency matrix A.\n",
    "\n",
    "    A[i, j] == 1 iff (i+1, j+1) is an edge (vertices in the file are 1-indexed).\n",
    "    Lines: 'c ...' comments, 'p edge <n> <m>' header, 'e <u> <v>' edges.\n",
    "    \"\"\"\n",
    "    n = None\n",
    "    edges = []\n",
    "    with open(file_path, 'r') as f:\n",
    "        for line in f:\n",
    "            parts = line.split()\n",
    "            if not parts:\n",
    "                continue\n",
    "            tag = parts[0]\n",
    "            if tag == 'c':\n",
    "                continue\n",
    "            elif tag == 'p':\n",
    "                # p edge n m   (some files use 'p col n m' / 'p clq n m')\n",
    "                n = int(parts[2])\n",
    "                m = int(parts[3])\n",
    "            elif tag == 'e':\n",
    "                edges.append((int(parts[1]), int(parts[2])))\n",
    "\n",
    "    if n is None:\n",
    "        raise ValueError(f\"no 'p' header found in {file_path}\")\n",
    "\n",
    "    A = np.zeros((n, n), dtype=np.uint8)\n",
    "    for u, v in edges:\n",
    "        if u == v:\n",
    "            continue  # ignore self-loops\n",
    "        A[u - 1, v - 1] = 1\n",
    "        A[v - 1, u - 1] = 1  # undirected -> symmetric\n",
    "\n",
    "    return A, n, m\n",
    "\n",
    "def build_clique_hamiltonian(A):\n",
    "    A = np.asarray(A, dtype=float)\n",
    "    return -A"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bcda93cd",
   "metadata": {},
   "source": [
    "### Functions to extract clique number from energy and statevector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58a60d5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def clique_number_from_energy(energy, sum_constraint=1.0):\n",
    "    R = float(sum_constraint)\n",
    "    return 1.0 / (1.0 + energy / (R * R))\n",
    "\n",
    "def extract_clique(x, A, tol=1e-6, scale=None):\n",
    "    \"\"\"Read a vertex set off the support of x and check it really is a clique.\n",
    "\n",
    "    `tol` is RELATIVE: v is in the support when x[v] > tol * scale, with scale\n",
    "    defaulting to max(x). \n",
    "    \n",
    "    Returns (vertices, is_clique). Vertices are 0-indexed into A; add 1 to get\n",
    "    the labels used in the DIMACS file.\n",
    "    \"\"\"\n",
    "    x = np.asarray(x, dtype=float)\n",
    "    if scale is None:\n",
    "        scale = x.max() if x.size else 0.0\n",
    "    if scale <= 0.0:\n",
    "        return np.array([], dtype=int), False\n",
    "    vertices = np.flatnonzero(x > tol * scale)\n",
    "    k = len(vertices)\n",
    "    if k == 0:\n",
    "        return vertices, False\n",
    "    sub = np.asarray(A, dtype=bool)[np.ix_(vertices, vertices)]\n",
    "    # A clique on k vertices has k*(k-1) ones in its (symmetric, zero-diagonal) block.\n",
    "    is_clique = bool(sub.sum() == k * (k - 1))\n",
    "    return vertices, is_clique\n",
    "\n",
    "def _extend_to_maximal(chosen, cand, Ab, x):\n",
    "    \"\"\"Grow `chosen` while candidates remain, highest weight first.\n",
    "    \"\"\"\n",
    "    while True:\n",
    "        idx = np.flatnonzero(cand)\n",
    "        if idx.size == 0:\n",
    "            return chosen\n",
    "        w = x[idx]\n",
    "        best = w.max()\n",
    "        ties = idx[w >= best - 1e-12 * max(1.0, abs(best))]\n",
    "        if ties.size > 1:\n",
    "            deg = Ab[np.ix_(ties, idx)].sum(1)\n",
    "            v = int(ties[np.argmax(deg)])\n",
    "        else:\n",
    "            v = int(ties[0])\n",
    "        chosen.append(v)\n",
    "        cand &= Ab[v]\n",
    "        cand[v] = False\n",
    "\n",
    "def greedy_clique_from_weights(x, A):\n",
    "    \"\"\"Repair step: greedily grow a genuine clique, taking vertices in order of x.\n",
    "    \"\"\"\n",
    "    Ab = np.asarray(A, dtype=bool)\n",
    "    x = np.asarray(x, dtype=float)\n",
    "    cand = np.ones(Ab.shape[0], dtype=bool)\n",
    "    return np.array(sorted(_extend_to_maximal([], cand, Ab, x)), dtype=int)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9cb2b5f",
   "metadata": {},
   "source": [
    "### Set up dirac-3S solver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5380ecc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def direct_dirac3(indices, coefficients, relaxation_schedule, sum_constraint, num_samples, ip_address):\n",
    "    eqc_client = EqcClient(ip_address=ip_address) #S9 \"172.18.41.45\", #S3 \"172.18.41.173\"\n",
    "    lock_id, start_ts, end_ts=eqc_client.wait_for_lock()\n",
    "    try:\n",
    "        result_dict = eqc_client.solve_sum_constrained(\n",
    "            lock_id=lock_id,\n",
    "            poly_indices = indices,\n",
    "            poly_coefficients = coefficients,\n",
    "            relaxation_schedule = relaxation_schedule,\n",
    "            sum_constraint = sum_constraint,\n",
    "            num_samples = num_samples,\n",
    "        )\n",
    "        total_time = np.array(result_dict[\"preprocessing_time\"])+np.array(result_dict[\"postprocessing_time\"])+np.array(result_dict[\"runtime\"])\n",
    "        print(f\"Total execution time(s):{total_time}\")\n",
    "    finally:\n",
    "        # release lock when finished using the device\n",
    "        lock_release_out = eqc_client.release_lock(lock_id=lock_id)\n",
    "    return result_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "824da4a4",
   "metadata": {},
   "source": [
    "### load instance file and solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ea5e9814",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "file loaded sucessfully\n"
     ]
    }
   ],
   "source": [
    "instance_dir = \"Instances/\"\n",
    "instance_name = \"keller4.clq\"\n",
    "instance_path = os.path.join(instance_dir, instance_name)\n",
    "try:\n",
    "    A, n, m = read_dimacs_clq(instance_path)\n",
    "    print(\"file loaded sucessfully\")\n",
    "except FileNotFoundError:\n",
    "    print(f\"{instance_path} does not exist\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e70d4ef2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "keller4: n = 171 vertices, m = 9435 edges (declared)\n",
      "A shape: (171, 171)\n",
      "edges in A: 9435\n",
      "symmetric: True\n",
      "self-loops: 0\n",
      "density: 0.64912\n"
     ]
    }
   ],
   "source": [
    "name = os.path.splitext(instance_name)[0]\n",
    "print(f\"{name}: n = {n} vertices, m = {m} edges (declared)\")\n",
    "print(f\"A shape: {A.shape}\")\n",
    "print(f\"edges in A: {int(A.sum()) // 2}\")\n",
    "print(f\"symmetric: {np.array_equal(A, A.T)}\")\n",
    "print(f\"self-loops: {int(np.trace(A))}\")\n",
    "print(f\"density: {A.sum() / (n * (n - 1)):.5f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0afa00e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "c = np.zeros(n)\n",
    "H = build_clique_hamiltonian(A)\n",
    "\n",
    "sum_constraint = 1\n",
    "num_samples = 100\n",
    "relaxation_schedule = 1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e97bf1c0",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:root:Max precision for EQC device is float32 input type was dtype float64. Input matrix will be rounded\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Submitting job to Dirac-3S\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/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\n",
      "  warnings.warn(warn_dtype_msg, Warning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total execution time(s):[2.62123575 2.58713802 2.01370343 2.7587886  2.04803766 2.67565884\n",
      " 2.23382308 2.79223963 2.67314126 2.64704813 2.84181871 2.49983934\n",
      " 2.79587569 2.90989884 2.82071676 2.66952636 2.52615025 1.86245462\n",
      " 2.70607432 2.54609944 2.66275682 2.76443251 2.69481905 2.16523417\n",
      " 2.80123534 2.7616445  2.03196877 2.04042424 2.80532928 2.7974691\n",
      " 2.62519749 2.03304568 2.78892185 2.24321626 2.80168899 2.83994545\n",
      " 2.62951382 2.03748437 2.54187473 2.63020476 2.86996313 2.68325921\n",
      " 2.79529806 2.86349048 1.97953007 2.16283818 2.18016803 2.97041192\n",
      " 2.82335903 2.07952028 2.68332769 2.79461773 2.54081376 2.40717832\n",
      " 2.45432718 2.01539499 2.68938633 2.20637588 2.57783373 2.75900568\n",
      " 2.50343429 2.69706388 2.68904926 2.06637112 2.11271995 2.66742732\n",
      " 2.26051008 2.03095476 2.69906934 2.74248256 2.79104366 1.99890017\n",
      " 2.78397586 2.7336476  2.65583034 2.90891715 2.6113561  2.27965211\n",
      " 2.732641   2.80390247 2.15183202 2.81881728 2.88916387 2.75707759\n",
      " 2.90165343 2.16630219 2.73159422 2.72164734 2.79233469 2.78964468\n",
      " 2.03095688 1.87199597 2.63561304 2.66995545 2.76049509 2.82846047\n",
      " 2.76028766 2.23443906 2.57353481 2.42797192]\n",
      "Dirac-3 Run complete. Response received.\n",
      "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]\n",
      "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]\n"
     ]
    }
   ],
   "source": [
    "poly_indices, poly_coefficients = convert_hamiltonian_to_poly_format(\n",
    "                    linear_terms=c,\n",
    "                    quadratic_terms=H,   # the util mutates its input\n",
    "                )\n",
    "\n",
    "dirac_start = time.time()\n",
    "print(\"Submitting job to Dirac-3S\")\n",
    "response = direct_dirac3(indices=poly_indices,\n",
    "                        coefficients=poly_coefficients,\n",
    "                        relaxation_schedule=relaxation_schedule,\n",
    "                        sum_constraint=sum_constraint,\n",
    "                        ip_address=ip_address,\n",
    "                        num_samples=num_samples)\n",
    "print(f\"Dirac-3 Run complete. Response received.\")\n",
    "energies = response['energy']\n",
    "solutions = response['solution']\n",
    "runtimes = response['runtime']\n",
    "preprocessing_time = response['preprocessing_time']\n",
    "postprocessing_time = response['postprocessing_time']\n",
    "\n",
    "print(f\"Dirac-3 all samples energies:{[round(e, 2) for e in energies]}\")\n",
    "print(f\"Dirac-3 all samples runtimes:{[round(t, 2) for t in runtimes]}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "bc966147",
   "metadata": {},
   "outputs": [],
   "source": [
    "best = min(zip(energies, solutions, runtimes), key=lambda x: x[0])\n",
    "dirac_energy, best_solution, best_runtime = best\n",
    "omega = clique_number_from_energy(dirac_energy, sum_constraint)\n",
    "support, isclique = extract_clique(best_solution, A)\n",
    "\n",
    "# Round every restart, not just the lowest-energy one.\n",
    "cliques = [greedy_clique_from_weights(sol, A) for sol in solutions]\n",
    "greedy = max(cliques, key=len)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "cf2632c6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "number of samples:100, relaxatrion schedule:1\n",
      "best energy over 100 samples:-0.888889\n",
      "support size:15, support is a clique:False\n",
      "clique vertices (1-indexed):[6, 22, 30, 39, 49, 74, 79, 105, 112]\n",
      "          energy   omega_est   clique    time (s)\n",
      "       -0.888889       9.000        9       1.986\n"
     ]
    }
   ],
   "source": [
    "print(f\"number of samples:{num_samples}, relaxatrion schedule:{relaxation_schedule}\")\n",
    "print(f\"best energy over {num_samples} samples:{dirac_energy:.6f}\")\n",
    "print(f\"support size:{len(support)}, support is a clique:{isclique}\")\n",
    "print(f\"clique vertices (1-indexed):{(greedy + 1).tolist()}\")\n",
    "print(f\"{'energy':>16}{'omega_est':>12}{'clique':>9}{'time (s)':>12}\")\n",
    "print(f\"{dirac_energy:>16.6f}{omega:>12.3f}{len(greedy):>9}{best_runtime:>12.3f}\")\n",
    "\n",
    "# Never report a set that is not actually a clique.\n",
    "assert extract_clique(np.isin(np.arange(n), greedy).astype(float), A)[1], \\\n",
    "    \"greedy repair returned a non-clique\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97982403",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}