Max_Cut.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "positive-possibility",
   "metadata": {},
   "source": [
    "# MAX-CUT for graphs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "declared-curve",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import math\n",
    "import pulp\n",
    "from random import *\n",
    "from pyvis.network import Network as netviz\n",
    "import networkx as net"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adjustable-incidence",
   "metadata": {},
   "source": [
    "### Generate a G(n,p)-random graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "id": "sharing-improvement",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100\n",
      "228\n"
     ]
    }
   ],
   "source": [
    "N,p = 100, 0.05\n",
    "G = net.gnp_random_graph(N,p)\n",
    "print(len(G.nodes))\n",
    "print(len(G.edges))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "broken-dining",
   "metadata": {},
   "source": [
    "## Create the adjacency matrix of G\n",
    "The function `adjacency_matrix()` returns a sparse matrix (i.e. only the coordinates of non-null entries). We turn it into a dense matrix before we can use it with the function `todense()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "dedicated-garlic",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = net.adjacency_matrix(G).todense()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "offshore-color",
   "metadata": {},
   "source": [
    "## Define a function to calculate the cut number for given `x`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "id": "about-twenty",
   "metadata": {},
   "outputs": [],
   "source": [
    "def CUT(G,x):\n",
    "    A = net.adjacency_matrix(G).todense()\n",
    "    n = len(G.nodes)\n",
    "    cut = 0\n",
    "    if len(x) != n:\n",
    "        print(\"ERROR length of x does not match dimension of G\")\n",
    "        return(cut)\n",
    "    else:\n",
    "        for i in range(1,n):\n",
    "            for j in range(1,n):\n",
    "                cut += 1/4*A.A[i][j]*(1-x[i]*x[j])\n",
    "        return(cut)  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "lasting-conclusion",
   "metadata": {},
   "source": [
    "## Generate $X \\sim \\mathcal{U} (\\{-1,1\\})$ and calculate CUT(G,X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "id": "irish-mitchell",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "121.0"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def rand_part(G):\n",
    "    \"return a random partition of V(G) represented by a vector with independent symmetric Bernoulli entries\"\n",
    "    X = []\n",
    "    for i in range(N):\n",
    "        X.append(1-2*round(random()))\n",
    "    return(X)\n",
    "X = rand_part(G)\n",
    "CUT(G,X)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "lesbian-oliver",
   "metadata": {},
   "source": [
    "## Define a (brute-force) function to calculate MAXCUT\n",
    "*WARNING* the runtime is $O(2^n)$ and _will_ take long for Graph sizes above $N=15$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "id": "operational-confidentiality",
   "metadata": {},
   "outputs": [],
   "source": [
    "def iterate(x):\n",
    "    x[0] = x[0]*(-1)\n",
    "    for i in range(len(x)):\n",
    "        if (x[i] == 1 ):\n",
    "            try:\n",
    "                x[i+1] = x[i+1]*(-1)\n",
    "            except IndexError:\n",
    "                break\n",
    "        else:\n",
    "            break\n",
    "    return(x)\n",
    "    \n",
    "    \n",
    "def MAXCUT(G):\n",
    "    n = len(G.nodes)\n",
    "    x = [1 for i in range(n)]\n",
    "    x0 = x\n",
    "    cut = 0\n",
    "    h = 0\n",
    "    for i in range (2**n):\n",
    "        h = CUT(G,x)\n",
    "        if h > cut:\n",
    "            cut = h\n",
    "            x0 = x\n",
    "            print(f'found a better CUT. NEW CUT# = {cut}')\n",
    "        x = iterate(x)\n",
    "    return(cut,x0)\n",
    "        \n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "incomplete-level",
   "metadata": {},
   "source": [
    "## A $(0.5-\\varepsilon)$-approximation algorithm for MAXCUT\n",
    "Let $\\varepsilon>0$ we start with a random partition of $V$, represented by a vector $X$ with independent symmetric Bernoulli entries. Then flip the entries of $X$ one by one and check if the cut number of the new cutset `CUT(G,X)` has improved. If so, update $X$ and continue until `CUT(G,X)` $\\geq (0.5-\\varepsilon)|E|$. As long as there are vertices where less than half of the incident edges are cut, the cut number is improved by flipping such a vertex. Hence, the algorithm improves can find a vertex to flip as long as `CUT(G,X)` $< 0.5 |E|$. Furthermore, the cut number increases by at least 1 with every flip, so the the number of flips is upper bounded by $|E|$.\n",
    "\n",
    "**Exercise** Implement this algorithm!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "id": "assigned-coordinator",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0, 1, 2, 3, 4, 6, 7, 9, 11, 12, 15, 16, 18, 19, 21, 22, 26, 27, 29, 31, 32, 33, 34, 35, 36, 38, 40, 42, 46, 47, 48, 52, 54, 55, 56, 58, 62, 63, 67, 69, 71, 72, 73, 77, 78, 79, 80, 81, 87, 89, 91, 96, 99]\n",
      "164.0\n"
     ]
    }
   ],
   "source": [
    "# Stefans solution\n",
    "\n",
    "def find_lower_edges(G,v):\n",
    "    \"\"\"helper function for maxcut_approx_Stefan\"\"\"\n",
    "    r = []\n",
    "    for u in G.neighbors(v):\n",
    "        if u < v:\n",
    "            r.append(u)\n",
    "    return(r)\n",
    "\n",
    "def v1_or_v2(G,V1,v):\n",
    "    \"\"\"helper function for maxcut_approx_Stefan\"\"\"\n",
    "    ln = find_lower_edges(G,v)\n",
    "    cnt = 0\n",
    "    r = False\n",
    "    for u in ln:\n",
    "        if u in V1:\n",
    "            cnt += 1\n",
    "    if cnt <= len(ln)/2:\n",
    "        r = True\n",
    "    return(r)\n",
    "\n",
    "def maxcut_approx_Stefan(G):\n",
    "    N = len(G.nodes())\n",
    "    V1 = [0]\n",
    "    for v in range(1,N):\n",
    "        if v1_or_v2(G,V1,v):\n",
    "            V1.append(v)\n",
    "    X = [-1 for i in range(N)]\n",
    "    print(V1)\n",
    "    for v in V1:\n",
    "        X[v] = 1\n",
    "    return(X)\n",
    "    \n",
    "X = maxcut_approx_Stefan(G)\n",
    "\n",
    "print(CUT(G,X))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "id": "076baf0d-1d6c-4bce-8232-bb7823dac818",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "175.0"
      ]
     },
     "execution_count": 93,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Thomas' solution\n",
    "def maxcut_approx_Thomas(G):\n",
    "    X = rand_part(G)\n",
    "    Cutnumber = CUT(G, X)\n",
    "    counter = 0 \n",
    "    while counter < N:\n",
    "        for i in range(N):\n",
    "            X[i] = -X[i]\n",
    "            if CUT(G, X) > Cutnumber:\n",
    "                Cutnumber = CUT(G, X)\n",
    "                counter = 0 \n",
    "            else:\n",
    "                X[i] = -X[i]\n",
    "                counter = counter + 1\n",
    "    return(X)\n",
    "\n",
    "CUT(G, maxcut_approx_Thomas(G)) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "finished-forward",
   "metadata": {},
   "source": [
    "## Draw the partitioned Graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "id": "australian-inside",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"800\"\n",
       "            height=\"800\"\n",
       "            src=\"net.html\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.IFrame at 0x7f18778448d0>"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def draw_part_graph(G,x):\n",
    "    n = len(G.nodes)\n",
    "    if len(x) != n:\n",
    "        print(\"ERROR length of x does not match dimension of G\")\n",
    "        return\n",
    "    for v in G.nodes():\n",
    "        if x[v] == 1:\n",
    "            G.nodes[v][\"color\"] = \"red\"\n",
    "        else:\n",
    "            G.nodes[v][\"color\"] = \"blue\"\n",
    "    edges = [e for e in G.edges]\n",
    "    for e in edges:\n",
    "        if x[e[0]]*x[e[1]] == -1:\n",
    "            G.edges[e][\"color\"] = \"red\"\n",
    "        else:\n",
    "            G.edges[e][\"color\"] = \"black\"\n",
    "    viz = netviz(height=800, width=800, notebook=True)\n",
    "    #viz.barnes_hut()\n",
    "    #viz.repulsion()\n",
    "    #viz.hrepulsion()\n",
    "    viz.force_atlas_2based()\n",
    "    viz.from_nx(G)\n",
    "    viz.show(\"net.html\")    \n",
    "    return(viz)\n",
    "\n",
    "viz = draw_part_graph(G,X)\n",
    "viz.show(\"net.html\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "accredited-exploration",
   "metadata": {},
   "source": [
    "## Calculate the MAXCUT-norm of G"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "nonprofit-parker",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "found a better CUT. NEW CUT# = 6.0\n",
      "found a better CUT. NEW CUT# = 8.0\n",
      "found a better CUT. NEW CUT# = 11.0\n",
      "found a better CUT. NEW CUT# = 13.0\n",
      "found a better CUT. NEW CUT# = 14.0\n",
      "found a better CUT. NEW CUT# = 17.0\n",
      "found a better CUT. NEW CUT# = 19.0\n",
      "found a better CUT. NEW CUT# = 22.0\n",
      "found a better CUT. NEW CUT# = 24.0\n",
      "found a better CUT. NEW CUT# = 27.0\n",
      "found a better CUT. NEW CUT# = 29.0\n",
      "found a better CUT. NEW CUT# = 30.0\n",
      "found a better CUT. NEW CUT# = 32.0\n",
      "found a better CUT. NEW CUT# = 35.0\n",
      "found a better CUT. NEW CUT# = 37.0\n",
      "found a better CUT. NEW CUT# = 40.0\n",
      "found a better CUT. NEW CUT# = 42.0\n",
      "found a better CUT. NEW CUT# = 44.0\n",
      "found a better CUT. NEW CUT# = 47.0\n",
      "found a better CUT. NEW CUT# = 49.0\n",
      "found a better CUT. NEW CUT# = 50.0\n",
      "found a better CUT. NEW CUT# = 52.0\n",
      "found a better CUT. NEW CUT# = 54.0\n",
      "found a better CUT. NEW CUT# = 55.0\n",
      "found a better CUT. NEW CUT# = 57.0\n",
      "found a better CUT. NEW CUT# = 59.0\n",
      "found a better CUT. NEW CUT# = 60.0\n",
      "found a better CUT. NEW CUT# = 62.0\n",
      "found a better CUT. NEW CUT# = 63.0\n",
      "found a better CUT. NEW CUT# = 64.0\n",
      "found a better CUT. NEW CUT# = 65.0\n",
      "found a better CUT. NEW CUT# = 66.0\n",
      "found a better CUT. NEW CUT# = 67.0\n",
      "found a better CUT. NEW CUT# = 68.0\n",
      "found a better CUT. NEW CUT# = 69.0\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-12-a7932e71db64>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mcut\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx0\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mMAXCUT\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mG\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcut\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-9-61d191b0327f>\u001b[0m in \u001b[0;36mMAXCUT\u001b[0;34m(G)\u001b[0m\n\u001b[1;32m     19\u001b[0m     \u001b[0mh\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     20\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 21\u001b[0;31m         \u001b[0mh\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mCUT\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mG\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     22\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mh\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0mcut\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     23\u001b[0m             \u001b[0mcut\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mh\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-7-f01e4a622496>\u001b[0m in \u001b[0;36mCUT\u001b[0;34m(G, x)\u001b[0m\n\u001b[1;32m      9\u001b[0m         \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     10\u001b[0m             \u001b[0;32mfor\u001b[0m \u001b[0mj\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 11\u001b[0;31m                 \u001b[0mcut\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m4\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mA\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mA\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mj\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mj\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     12\u001b[0m         \u001b[0;32mreturn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcut\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/tljh/user/lib/python3.7/site-packages/numpy/matrixlib/defmatrix.py\u001b[0m in \u001b[0;36mA\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    863\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    864\u001b[0m         \"\"\"\n\u001b[0;32m--> 865\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__array__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    866\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    867\u001b[0m     \u001b[0;34m@\u001b[0m\u001b[0mproperty\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "cut, x0 = MAXCUT(G)\n",
    "print(cut)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f20b16f-0878-4f76-9cb1-f6c23a96c7dc",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}