Graph_Isomorphism_Lecture.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Graph Isomorphism Problem\n",
    "\n",
    "The Graph Isomorphism Problem asks whether two given finite graphs are isomorphic.\n",
    "An isomorphism of graphs $G_1$ and $G_2$ is a bijection (one-to-one mapping) between the vertex sets of $G_1$ and $G_2$,\n",
    "such that the adjacency structure is preserved.\n",
    "\n",
    "The Graph Isomorphism Problem and variations of it have many applications in areas such as:\n",
    "- image processing\n",
    "- protein structure\n",
    "- computer systems\n",
    "- malware detection\n",
    "- chemical bond structure\n",
    "- social networks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<table class=\"wikitable\" style=\"margin: 1em auto 1em auto\">\n",
    "<tbody><tr>\n",
    "<th>Graph G1\n",
    "</th>\n",
    "<th>Graph G2\n",
    "</th>\n",
    "<th>An isomorphism<br>between G1 and G2\n",
    "</th></tr>\n",
    "<tr>\n",
    "<td style=\"padding-left:2em;padding-right:2em;\"><a href=\"/wiki/File:Graph_isomorphism_a.svg\" class=\"image\"><img alt=\"Graph isomorphism a.svg\" src=\"//upload.wikimedia.org/wikipedia/commons/thumb/9/9a/Graph_isomorphism_a.svg/100px-Graph_isomorphism_a.svg.png\" decoding=\"async\" width=\"100\" height=\"200\" srcset=\"//upload.wikimedia.org/wikipedia/commons/thumb/9/9a/Graph_isomorphism_a.svg/150px-Graph_isomorphism_a.svg.png 1.5x, //upload.wikimedia.org/wikipedia/commons/thumb/9/9a/Graph_isomorphism_a.svg/200px-Graph_isomorphism_a.svg.png 2x\" data-file-width=\"200\" data-file-height=\"400\"></a>\n",
    "</td>\n",
    "<td style=\"padding-left:1em;padding-right:1em;\"><a href=\"/wiki/File:Graph_isomorphism_b.svg\" class=\"image\"><img alt=\"Graph isomorphism b.svg\" src=\"//upload.wikimedia.org/wikipedia/commons/thumb/8/84/Graph_isomorphism_b.svg/210px-Graph_isomorphism_b.svg.png\" decoding=\"async\" width=\"210\" height=\"210\" srcset=\"//upload.wikimedia.org/wikipedia/commons/thumb/8/84/Graph_isomorphism_b.svg/315px-Graph_isomorphism_b.svg.png 1.5x, //upload.wikimedia.org/wikipedia/commons/thumb/8/84/Graph_isomorphism_b.svg/420px-Graph_isomorphism_b.svg.png 2x\" data-file-width=\"400\" data-file-height=\"400\"></a>\n",
    "</td>\n",
    "    <td align=\"center\" style=\"background-color:white;\">\n",
    "\n",
    "<div style=\"padding:1em\"><i>f</i>(<i>a</i>) = 1\n",
    "</div><div style=\"padding:0.5em\"><i>f</i>(<i>b</i>) = 6\n",
    "</div><div style=\"padding:0.5em\"><i>f</i>(<i>c</i>) = 8\n",
    "</div><div style=\"padding:0.5em\"><i>f</i>(<i>d</i>) = 3\n",
    "</div><div style=\"padding:0.5em\"><i>f</i>(<i>g</i>) = 5\n",
    "</div><div style=\"padding:0.5em\"><i>f</i>(<i>h</i>) = 2\n",
    "</div><div style=\"padding:0.5em\"><i>f</i>(<i>i</i>) = 4\n",
    "</div><div style=\"padding:0.5em\"><i>f</i>(<i>j</i>) = 7\n",
    "</div>\n",
    "</td></tr></tbody></table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem formulation - Objective and Constraints\n",
    "\n",
    "## Decision Variables\n",
    "Let's define a binary decision variable $x_{i,j}$ to determine if vertex $i$ of one graph ($G_1$) maps to vertex $j$ of another graph ($G_2$).\n",
    "\n",
    "If $x_{i,j} = 1$, then we know $f(i) = j$ and $f^{-1}(j) = i$.\n",
    "\n",
    "## Objective Function\n",
    "\n",
    "We want to define the objective function (or Hamiltonian $H$) as a series of penalties for pairs of vertices that are connected in one graph and not in the other.\n",
    "\n",
    "- Penalize all \"bad edge mappings\" (edge in $E_{G_1}$, but not in $E_{G_2}$, vice versa)\n",
    "\n",
    "$$\\Large H = \\sum_{ij \\notin E_{G_1}} \\sum_{uv \\in E_{G_2}} x_{i, u}x_{j,v}+\\sum_{ij \\in E_{G_1}} \\sum_{uv \\notin E_{G_2}} x_{i, u}x_{j,v}$$\n",
    "\n",
    "For large, sparse graphs, this can produce a large number of biases because we are penalizing missing edges.\n",
    "\n",
    "\n",
    "## Constraints\n",
    "\n",
    "Based on the definition of the decision variable $x_{i,j}$, each vertex $i$ in graph $G_1$ can only be mapped to one vertex in the second graph and vice versa.\n",
    "\n",
    "$$\\Large\\sum_{j \\in V_{G_2}} x_{i,j} = 1 , \\forall i \\in G_1$$\n",
    "\n",
    "$$\\Large\\sum_{i \\in V_{G_1}} x_{i,j} = 1 , \\forall j \\in G_2 $$\n",
    "\n",
    "\n",
    "## Discrete Quadratic Model (DQM)\n",
    "Discrete quadratic model naturally handles one-hot constraints. We can think of different values of $j \\in V_{G_2}$ as cases of each discrete variable $x_i$, such that if the binary representation $x_{i,j} = 1$, then the discrete representation is $x_i = j$. Conversely, $x_i = j$ implies not only that $x_{i,j}=1$ but also that all other $x_{i,k}=0$ if $k \\neq j$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Code to solve an instance of the Graph Isomorphism Problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import networkx as nx\n",
    "\n",
    "num_nodes = 10\n",
    "graph_seed = 42\n",
    "\n",
    "# Two graphs we know to be isomorphic by design\n",
    "g_1 = nx.gnp_random_graph(n=num_nodes, p=0.5, seed=graph_seed)\n",
    "g_1 = nx.relabel_nodes(g_1, {i: chr(i + ord('a')) for i in range(num_nodes)})\n",
    "\n",
    "g_2 = nx.gnp_random_graph(n=num_nodes, p=0.5, seed=graph_seed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, ax = plt.subplots(1, 2, figsize=(10, 6))\n",
    "nx.draw(g_1, ax=ax[0], with_labels=True, font_color='w', font_size=18, node_size=500)\n",
    "nx.draw(g_2, ax=ax[1], with_labels=True, font_color='w', font_size=18, node_size=500)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dimod import DQM\n",
    "\n",
    "dqm = DQM()\n",
    "\n",
    "for node in g_1.nodes:\n",
    "    dqm.add_variable(num_nodes, label=node)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Objective"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from itertools import combinations\n",
    "\n",
    "penalty = 1\n",
    "for (i, j) in g_1.edges:\n",
    "    for u, v in combinations(g_2.nodes, r=2):\n",
    "        if (u, v) not in g_2.edges:\n",
    "            dqm.set_quadratic_case(i, u, j, v, penalty)\n",
    "            dqm.set_quadratic_case(i, v, j, u, penalty)\n",
    "\n",
    "for (u, v) in g_2.edges:\n",
    "    for i, j in combinations(g_1.nodes, r=2):\n",
    "        if (i, j) not in g_1.edges:\n",
    "            dqm.set_quadratic_case(i, u, j, v, penalty)\n",
    "            dqm.set_quadratic_case(i, v, j, u, penalty)\n",
    "\n",
    "print('Total number of variables: ', dqm.num_variables())\n",
    "print('Number of cases per variable: ', dqm.num_cases()//dqm.num_variables())\n",
    "print('Total number of pairwise interactions: ', dqm.num_case_interactions())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Constraints\n",
    "\n",
    "Since the Discrete Quadratic Model naturally handles one-hot constraints, variables representing vertices in $G_1$ will only assume one case representing a vertex in $G_2$.\n",
    "\n",
    "We still need to ensure that the variables representing each vertex in $G_1$ map to unique vertices in $G_2$, which is our bijection constraint above.\n",
    "\n",
    "To encode this constraint, we must convert it into an objective: \n",
    "$$ \\Large (\\sum_{j \\in V_{G_2}} x_{i,j} - 1)^2 $$\n",
    "\n",
    "At this point we would expand the expression, however the Discrete Quadratic Model has a nice feature for these constraints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for u in g_2.nodes:\n",
    "    dqm.add_linear_equality_constraint(\n",
    "        [(i, u, 1.0) for i in g_1.nodes],\n",
    "        constant=-1,\n",
    "        lagrange_multiplier=penalty\n",
    "    )\n",
    "\n",
    "print('Total number of pairwise interactions: ', dqm.num_case_interactions())\n",
    "\n",
    "print('Offset: ', dqm.offset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve with Leap Hybrid DQM Sampler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dwave.system import LeapHybridDQMSampler\n",
    "\n",
    "sampler = LeapHybridDQMSampler()\n",
    "\n",
    "sampleset = sampler.sample_dqm(dqm, time_limit=5)\n",
    "sampleset.record.energy = dqm.energies(sampleset)\n",
    "\n",
    "bijection = sampleset.first.sample\n",
    "best_energy = sampleset.first.energy\n",
    "\n",
    "print('Bijection: ', bijection)\n",
    "print('Objective: ', best_energy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils import plot_graphs\n",
    "plot_graphs(g_1, g_2, bijection)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References and related demo\n",
    "\n",
    "- Intro imagery and information about problem from: https://en.wikipedia.org/wiki/Graph_isomorphism\n",
    "\n",
    "- D-Wave Examples circuit equivalence repository: https://github.com/dwave-examples/circuit-equivalence"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}