{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1a5c6312",
   "metadata": {},
   "source": [
    "# PGD Max Clique Benchmarking"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "daa813cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import os\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "999fb2d3",
   "metadata": {},
   "source": [
    "### Read data file and build objective functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3aa3a69e",
   "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": "1d3f375b",
   "metadata": {},
   "source": [
    "### Functions to extract clique number from energy and statevector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8064f8d0",
   "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": "7d3d2837",
   "metadata": {},
   "source": [
    "### Define multistart projected gradient descent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "345fbbb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def project_simplex(v, sum_constraint):\n",
    "    \"\"\"Exact Euclidean projection onto {x >= 0, sum(x) == R}.\n",
    "    \"\"\"\n",
    "    R = float(sum_constraint)\n",
    "    v = np.asarray(v, dtype=float)\n",
    "    u = np.sort(v)[::-1]\n",
    "    cumulative = np.cumsum(u) - R\n",
    "    ind = np.arange(1, v.size + 1)\n",
    "    rho = ind[u - cumulative / ind > 0][-1]\n",
    "    return np.maximum(v - cumulative[rho - 1] / rho, 0.0)\n",
    "\n",
    "\n",
    "def pgd(Q, c, sum_constraint, lr=0.01, max_iter=10**4, tol=1e-9, x0=None):\n",
    "    \"\"\"Projected gradient descent on {x >= 0, sum(x) == R}.\n",
    "    \"\"\"\n",
    "    Q = np.asarray(Q, dtype=float)\n",
    "    c = np.asarray(c, dtype=float).ravel()\n",
    "    n = Q.shape[0]\n",
    "    R = float(sum_constraint)\n",
    "\n",
    "    x = np.full(n, R / n) if x0 is None else np.asarray(x0, dtype=float)\n",
    "    x = project_simplex(x, R)\n",
    "    QT = Q + Q.T                                   # gradient of x'Qx is (Q+Q')x\n",
    "\n",
    "    it = 0\n",
    "    for it in range(max_iter):\n",
    "        grad = QT @ x + c\n",
    "        x_next = project_simplex(x - lr * grad, R)\n",
    "        delta = np.linalg.norm(x_next - x)\n",
    "        x = x_next                                 # keep the newer iterate\n",
    "        if delta < tol:\n",
    "            break\n",
    "\n",
    "    return x, float(x @ Q @ x + c @ x), it + 1\n",
    "\n",
    "\n",
    "def pgd_multistart(Q, c, sum_constraint, restarts=32, seed=0, **kwargs):\n",
    "    \"\"\"PGD from the simplex centre plus random starts.\n",
    "    \"\"\"\n",
    "    rng = np.random.default_rng(seed)\n",
    "    n = Q.shape[0]\n",
    "    solutions, energies, times, total_iters = [], [], [], 0\n",
    "    for r in range(restarts):\n",
    "        start = time.time()\n",
    "        x0 = None\n",
    "        if r > 0:\n",
    "            x0 = rng.random(n)\n",
    "            x0 *= sum_constraint / x0.sum()\n",
    "        x, energy, iters = pgd(Q, c, sum_constraint, x0=x0, **kwargs)\n",
    "        pgd_t = time.time()-start\n",
    "        solutions.append(x)\n",
    "        energies.append(energy)\n",
    "        times.append(pgd_t)\n",
    "        total_iters += iters\n",
    "    return solutions, energies, times, total_iters\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93c15915",
   "metadata": {},
   "source": [
    "### load instance file and solve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "48053939",
   "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": null,
   "id": "a70a3b8d",
   "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": 7,
   "id": "0fdc3aa6",
   "metadata": {},
   "outputs": [],
   "source": [
    "c = np.zeros(n)\n",
    "H = build_clique_hamiltonian(A)\n",
    "\n",
    "sum_constraint = 1\n",
    "restarts = 100\n",
    "learning_rate = 0.01\n",
    "max_iter = 2 * 10**4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "2860df89",
   "metadata": {},
   "outputs": [],
   "source": [
    "pgd_start = time.time()\n",
    "solutions, energies, times, total_iters = pgd_multistart(\n",
    "    Q=H,\n",
    "    c=c,\n",
    "    sum_constraint=sum_constraint,\n",
    "    restarts=restarts,\n",
    "    lr=learning_rate,\n",
    "    max_iter=max_iter,\n",
    ")\n",
    "pgd_time = time.time() - pgd_start"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f654846",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_energy = min(energies)\n",
    "best_solution = solutions[int(np.argmin(energies))]\n",
    "omega = clique_number_from_energy(best_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": null,
   "id": "84c8fd84",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "restarts:100, total PGD iterations:77765\n",
      "best energy over restarts:-0.888889\n",
      "support size:15, support is a clique:False\n",
      "clique vertices (1-indexed):[3, 24, 28, 37, 59, 77, 87, 105, 111]\n",
      "          energy   omega_est   clique    time (s)\n",
      "       -0.888889       9.000        9       0.735\n"
     ]
    }
   ],
   "source": [
    "\n",
    "print(f\"restarts:{restarts}, total PGD iterations:{total_iters}\")\n",
    "print(f\"best energy over restarts:{best_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\"{best_energy:>16.6f}{omega:>12.3f}{len(greedy):>9}{pgd_time:>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": "8a39dedb",
   "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
}