{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TODO\n",
    "- time-varying graphs\n",
    "- table of results\n",
    "- show combinatorics. why not NP hard. enumerate paths.\n",
    "- remove $\\mathcal{P}$. just describe a path with $x$ from the start?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Topic references\n",
    "- [Convex Optimization book](http://web.stanford.edu/~boyd/cvxbook/)\n",
    "- [Additional Exercise 15.4: Allocation of interdiction effort](http://web.stanford.edu/~boyd/cvxbook/bv_cvxbook_extra_exercises.pdf)\n",
    "- [Shortest path problem in graph theory](http://en.wikipedia.org/wiki/Shortest_path_problem)\n",
    "- [Dijkstra's_algorithm](http://en.wikipedia.org/wiki/Dijkstra%27s_algorithm)\n",
    "- [Iterated weighted $\\ell_1$ heuristic](http://stanford.edu/class/ee364b/lectures/l1_ext_slides.pdf)\n",
    "\n",
    "# Allocating interdiction effort to catch a smuggler\n",
    "\n",
    "We'll explore a game between a smuggler moving through a graph and a security team trying to catch him.\n",
    "\n",
    "The smuggler moves along a directed graph (cycles are allowed) with $n$ nodes and $m$ edges from node $0$ (the source) to node $n-1$ (the destination).\n",
    "Edge $j$ has evasion probability $p_j$, which is the probability that the smuggler moves along that edge undetected.\n",
    "The detection events on the edges are independent, so the probability that the smuggler makes it to the destination undetected is $\\prod_{j \\in \\mathcal{P}} p_j$, where $\\mathcal{P} \\subset \\lbrace 0,\\ldots, m-1\\rbrace$ is the set of edges on the smuggler's path.\n",
    "We assume that the smuggler knows the evasion probabilities and will take a path that maximizes the probability of making it to the destination node undetected.\n",
    "\n",
    "The security team is trying to catch the smuggler as he moves through the graph.\n",
    "They can control the evasion probabilities on the graph edges through, say, allocation of resources to the edges, assigning guards, etc.\n",
    "We'll keep this abstract for now, and just consider that security can modify the evasion probabilities subject to some constraints.\n",
    "Their goal is to choose evasion probabilities to minimize the chance that the smuggler gets through undetected.\n",
    "\n",
    "# Transformations and mathematical formulation\n",
    "\n",
    "We'll take the logarithm of all probabilities to make the problem amenable to convex optimization and simplify the remaining discussion.\n",
    "This transforms products of variables (generally not convex) into sums of variables (convex!).\n",
    "\n",
    "Let $c_j = \\log(p_j)$.\n",
    "As this is a one-to-one correspondence, we may refer to the 'edge detection probabilities' when we really mean to refer to their logarithm, $c$.\n",
    "This should be clear from context and hopefully won't cause any confusion.\n",
    "\n",
    "Given $c$ and a path $\\mathcal{P}$, let $P(c, \\mathcal{P})$ be the logarithm of the evasion probability of that path. That is, \n",
    "\n",
    "$$\n",
    "P(c, \\mathcal{P}) = \\log \\left( \\prod_{j \\in \\mathcal{P}} p_j \\right) = \\sum_{j \\in \\mathcal{P}} \\log(p_j) = \\sum_{j \\in \\mathcal{P}} c_j.\n",
    "$$\n",
    "\n",
    "Define $P(c, \\mathcal{P}) = -\\infty$ if $\\mathcal{P}$ is not a valid path from source to sink. We do this to implicitly encode the constraint that the smuggler's path must be valid.\n",
    "\n",
    "# Min/max game\n",
    "\n",
    "We can now interpret this problem as a min/max game. The smuggler is trying to maximize his chance of getting through undetected by choosing $\\mathcal{P}$, and security is trying to minimize the smuggler's chance by choosing $c$. That is, the smuggler is trying to maximize $P(c, \\mathcal{P})$ while the security team is trying to minimize it.\n",
    "\n",
    "We end up with a game with two 'moves'.\n",
    "Security 'moves' first by choosing the edge evasion probabilities via $c$.\n",
    "The smuggler 'moves' second (after $c$ is fixed) by choosing his path $\\mathcal{P}$.\n",
    "\n",
    "For security to make an optimal choice of edge evasion probabilities, they'll need to model what the smuggler will do when faced with any set of evasion probabilities. Thus, we'll let $P^\\mathrm{opt}(c)$ be the optimal value of the problem\n",
    "\n",
    "\\\\[\n",
    "\\begin{array}{ll}\n",
    "  \\mbox{maximize} & P(c, \\mathcal{P}) \\\\\n",
    "  \\mbox{subject to} & \\mathcal{P} \\mbox{ is a valid path}.\n",
    "\\end{array}\n",
    "\\\\]\n",
    "\n",
    "$P^\\mathrm{opt}(c)$ is the (logarithm of the) evasion probability of the smuggler if he picks an optimal path, knowing the values of $c$.\n",
    "\n",
    "The security team's goal is then to minimize this quantity, subject to some constraints on $c$, which we'll denote by the set $\\mathcal{C}$. Let $P^{\\star}$ be the optimal value of \n",
    "\n",
    "\\\\[\n",
    "\\begin{array}{ll}\n",
    "  \\mbox{minimize} & P^{\\mathrm{opt}}(c) \\\\\n",
    "  \\mbox{subject to} & c \\in \\mathcal{C}.\n",
    "\\end{array}\n",
    "\\\\]\n",
    "\n",
    "In the remainder, we'll first consider the smuggler's objective so as to obtain a convex formulation of \n",
    "$P^\\mathrm{opt}(c)$, which we can then use to analyze and optimize the security team's objective."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Helper functions\n",
    "\n",
    "The next cell contains some helper functions for generating the examples and plotting the results. This is the majority of the code. The actual CVXPY optimization code is only a few lines. We'll go over the optimization code later in the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%config InlineBackend.figure_format = 'pdf'\n",
    "#%config InlineBackend.figure_format = 'svg'\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import networkx as nx\n",
    "import cvxpy as cp\n",
    "\n",
    "from matplotlib import rcParams\n",
    "rcParams.update({'font.size': 16})\n",
    "\n",
    "def formGraph(N,mu,eta,buildings=None,seed=None,dummyNodes=False):\n",
    "    ''' Form N by N grid of nodes, perturb by mu and connect nodes within eta.\n",
    "        mu and eta are relative to 1/(N-1)\n",
    "        \n",
    "        buildings is a list of tuples (x,y,w,h)\n",
    "    '''\n",
    "    if seed is not None:\n",
    "        np.random.seed(seed)\n",
    "    \n",
    "    mu = mu/(N-1)\n",
    "    eta = eta/(N-1)\n",
    "    \n",
    "    # generate perterbed grid positions for the nodes\n",
    "    pos = [(i + mu*np.random.randn(), j + mu*np.random.randn())\\\n",
    "       for i in np.linspace(0,1,N) for j in np.linspace(1,0,N)]\n",
    "    \n",
    "    #select out nodes that end up inside buildings\n",
    "    if buildings is not None and buildings != []:\n",
    "        pos2 = []\n",
    "        for p in pos:\n",
    "            inside = False\n",
    "            for x,y,w,h in buildings:\n",
    "                if x <= p[0] and p[0] <= x + w and y <= p[1] and p[1] <= y + h:\n",
    "                    inside = True\n",
    "                    break\n",
    "            if not inside:\n",
    "                pos2.append(p)\n",
    "        pos = pos2\n",
    "        \n",
    "    # add dummy nodes for multiple entry/exit example\n",
    "    if dummyNodes:\n",
    "        pos = [(-.5,.5)] + pos + [(1.5,.5)]\n",
    "    \n",
    "    pos = dict(enumerate(pos))\n",
    "    n = len(pos)\n",
    "    \n",
    "    # connect nodes with edges\n",
    "    G = nx.random_geometric_graph(n,eta,pos=pos)\n",
    "    \n",
    "    # remove edges that cross buildings\n",
    "    if buildings is not None and buildings != []:\n",
    "        for e in list(G.edges()):\n",
    "            blocked = False\n",
    "            for x,y,w,h in buildings:\n",
    "                if intersect(pos[e[0]],pos[e[1]],x,x+w,y,y+h):\n",
    "                    blocked = True\n",
    "            if blocked:\n",
    "                G.remove_edge(*e)\n",
    "    \n",
    "    G = nx.DiGraph(G)\n",
    "    \n",
    "    # add edges connecting dummy nodes to nodes on left and right edges of grid\n",
    "    if dummyNodes:\n",
    "        for i in range(N):\n",
    "            G.add_edge(0,i+1)\n",
    "            G.add_edge(n-i-2,n-1)\n",
    "    \n",
    "    return G, pos\n",
    "\n",
    "def showPaths(G,pos,edgeProbs=1.0,path=None,visibleNodes=None,Gnodes=None,Rnodes=None,guards=None):\n",
    "    ''' Takes directed graph G, node positions pos, and edge probabilities.\n",
    "        Optionally uses path (a list of edge indices) to plot the smuggler's path.\n",
    "        \n",
    "        edgeProbd gives the probabilities for all the edges, including hidden ones.\n",
    "        \n",
    "        path includes all the edges, including the hidden ones\n",
    "        \n",
    "        Gnodes and Rnodes denote the source and destination nodes, to be plotted green\n",
    "        and red respectively.\n",
    "        \n",
    "        guards is a list of node indices for denoting guards with a black dot on the plot\n",
    "    '''\n",
    "    fig = plt.figure(figsize=(8,6))\n",
    "    ax = fig.add_subplot(111, aspect='equal')\n",
    "    \n",
    "    n = G.number_of_nodes()\n",
    "    if visibleNodes is None:\n",
    "        visibleNodes = G.nodes()\n",
    "        \n",
    "    # draw the regular interior nodes in the graph\n",
    "    nx.draw_networkx_nodes(G,pos,nodelist=visibleNodes, edgecolors='k',\n",
    "                           node_color='w',node_size=100,ax=ax)\n",
    "\n",
    "    # draw the origin and destination nodes\n",
    "    for nodes, color in zip([Gnodes,Rnodes],['g','r']):\n",
    "        for color2, alpha in zip(['w',color],[1,.2]):\n",
    "            nx.draw_networkx_nodes(G,pos,\n",
    "                           nodelist=nodes,\n",
    "                           node_color=color2,\n",
    "                           edgecolors='k',\n",
    "                           node_size=200,\n",
    "                           ax=ax,alpha=alpha)\n",
    "\n",
    "    # draw guard nodes\n",
    "    if guards is not None:\n",
    "        nx.draw_networkx_nodes(G,pos,nodelist=guards,node_color='k',\n",
    "                               node_size=100,ax=ax)\n",
    "      \n",
    "        \n",
    "    if path is None:\n",
    "        alpha = 1\n",
    "    else:\n",
    "        alpha = .15\n",
    "        \n",
    "    # start messing with edges\n",
    "    edge2ind = {e:i for i,e in enumerate(G.edges())}\n",
    "    ind2edge = {i:e for i,e in enumerate(G.edges())}\n",
    "    \n",
    "    # only display edges between non-dummy nodes\n",
    "    visibleEdges = [i for i in range(len(edge2ind)) if ind2edge[i][0] in visibleNodes and ind2edge[i][1] in visibleNodes]\n",
    "    \n",
    "    edgelist = [ind2edge[i] for i in visibleEdges]\n",
    "    \n",
    "    if isinstance(edgeProbs,float):\n",
    "        edgeProbs = [edgeProbs]*G.number_of_edges()\n",
    "        \n",
    "    p = [edgeProbs[i] for i in visibleEdges]\n",
    "\n",
    "    # draw edges of graph, make transparent if we're drawing a path over them\n",
    "    edges = nx.draw_networkx_edges(G,pos,edge_color=p,width=4,\n",
    "                                   edge_cmap=plt.cm.RdYlGn,arrows=False,edgelist=edgelist,edge_vmin=0.0,\n",
    "                                   edge_vmax=1.0,ax=ax,alpha=alpha)\n",
    "\n",
    "    # draw the path, only between visible nodes\n",
    "    if path is not None:\n",
    "        visiblePath = [i for i in path if ind2edge[i][0] in visibleNodes and ind2edge[i][1] in visibleNodes]\n",
    "        path_pairs = [ind2edge[i] for i in visiblePath]\n",
    "        path_colors = [edgeProbs[i] for i in visiblePath]\n",
    "        edges = nx.draw_networkx_edges(G,pos,edge_color=path_colors,width=8,\n",
    "                                       edge_cmap=plt.cm.RdYlGn,edgelist=path_pairs,arrows=False,edge_vmin=0.0,\n",
    "                                   edge_vmax=1.0)\n",
    "    \n",
    "    fig.colorbar(edges,label='Evasion probability')\n",
    "\n",
    "    ax.axis([-0.05,1.05,-0.05,1.05])\n",
    "    #ax.axis('tight')\n",
    "    #ax.axis('equal')\n",
    "    ax.axis('off')\n",
    "    \n",
    "    return fig, ax\n",
    "    \n",
    "def optPath(G,p):\n",
    "    ''' Find the optimal smuggler's path in directed graph G\n",
    "        with edge evasion probabilities p and dictionary\n",
    "        edge2ind bringing node pairs to edge indices\n",
    "    '''\n",
    "    edge2ind = {e:i for i,e in enumerate(G.edges())}\n",
    "    H = G.copy()\n",
    "    p = np.minimum(p,1)\n",
    "    w = -np.log(p)\n",
    "    n = H.number_of_nodes()\n",
    "    \n",
    "    for i in H:\n",
    "        for j in H[i]:\n",
    "            H[i][j]['w'] = w[edge2ind[(i,j)]]\n",
    "            \n",
    "    path = nx.shortest_path(H,0,n-1,'w')\n",
    "    #path = nx.astar_path(H,0,n-1,weight='w')\n",
    "            \n",
    "    foo = [edge2ind[(i,j)] for i,j in zip(path[:-1],path[1:])]\n",
    "    x = np.zeros_like(p)\n",
    "    x[foo] = 1\n",
    "    return x\n",
    "\n",
    "def intersect(p1,p2,xmin,xmax,ymin,ymax):\n",
    "    '''determine if a rectangle given by xy limits blocks the line of sight between p1 and p2'''\n",
    "\n",
    "    block = False\n",
    "    \n",
    "    # if either point inside block\n",
    "    for p in [p1,p1]:\n",
    "        if xmin <= p[0] and p[0] <= xmax and ymin <= p[1] and p[1] <= ymax:\n",
    "            return True\n",
    "    \n",
    "    # if the two points are equal at this stage, then they are outside the block\n",
    "    if p1[0] == p2[0] and p1[1] == p2[1]:\n",
    "        return False\n",
    "    \n",
    "    \n",
    "    if p2[0] != p1[0]:\n",
    "        for x in [xmin,xmax]:\n",
    "            alpha = (x-p1[0])/(p2[0] - p1[0])\n",
    "            y = p1[1] + alpha*(p2[1] - p1[1])\n",
    "\n",
    "            if 0 <= alpha and alpha <= 1 and ymin <= y and y <= ymax:\n",
    "                return True\n",
    "            \n",
    "    if p2[1] != p1[1]:\n",
    "        for y in [ymin,ymax]:\n",
    "            alpha = (y-p1[1])/(p2[1] - p1[1])\n",
    "            x = p1[0] + alpha*(p2[0] - p1[0])\n",
    "\n",
    "            if 0 <= alpha and alpha <= 1 and xmin <= x and x <= xmax:\n",
    "                return True\n",
    "        \n",
    "    return False\n",
    "\n",
    "def getGuardEffects(G,pos,guardIdxs,buildings=None,dummyNodes=None):\n",
    "    ''' for n guards and m edges, return an m by n matrix giving the\n",
    "        effect of each guard on the evasion probabilities of each edge\n",
    "        in the graph.\n",
    "        \n",
    "        Ignore dummy nodes, if any. Guards cannot see through buildings.\n",
    "        \n",
    "        guardIdxs is a list of the node indices where we would consider placing\n",
    "        guards\n",
    "        \n",
    "        Return the log of the detection probabilities for each guard.\n",
    "    \n",
    "    '''\n",
    "    def evadeProb(x,radius=.2):\n",
    "        r = np.linalg.norm(x)/radius\n",
    "        return min(r+.1,1)\n",
    "    \n",
    "    m = G.number_of_edges()\n",
    "    if buildings is None:\n",
    "        buildings = []\n",
    "    if dummyNodes is None:\n",
    "        dummyNodes = []\n",
    "    \n",
    "    \n",
    "    ind2edge = {i:e for i,e in enumerate(G.edges())}\n",
    "    edgeCenters = []\n",
    "    for e in G.edges():\n",
    "        edgeCenters.append((np.array(pos[e[0]]) + np.array(pos[e[1]]))/2)\n",
    "        \n",
    "    edgeCenters = np.array(edgeCenters)\n",
    "    numGuards = len(guardIdxs)\n",
    "    edgeVals = np.zeros((m,numGuards))\n",
    "    \n",
    "    for i,gid in enumerate(guardIdxs):\n",
    "        for j in range(m):\n",
    "            blocked = False\n",
    "            for x,y,w,h in buildings:\n",
    "                if intersect(edgeCenters[j,:],pos[gid],x,x+w,y,y+h):\n",
    "                    blocked = True\n",
    "                    break\n",
    "            e = ind2edge[j]\n",
    "            if e[0] in dummyNodes or e[1] in dummyNodes:\n",
    "                blocked = True\n",
    "            if not blocked:\n",
    "                edgeVals[j,i] += np.log(evadeProb(pos[gid]-edgeCenters[j,:],.3))\n",
    "    return edgeVals\n",
    "\n",
    "def get_a(G,seed=None):\n",
    "    '''\n",
    "    Generate a random value in [0,1] for each edge in the graph.\n",
    "    For directed graphs, directed edges between the same nodes should have the same value\n",
    "    '''\n",
    "    if seed is not None:\n",
    "        np.random.seed(seed)\n",
    "    m = G.number_of_edges()\n",
    "    a = np.zeros((m))\n",
    "    edge2ind = {e:i for i,e in enumerate(G.edges())}\n",
    "    for e in G.edges():\n",
    "        a[edge2ind[e]] = np.random.rand()\n",
    "        e2 = (e[1],e[0])\n",
    "        if e2 in edge2ind:\n",
    "            a[edge2ind[e2]] = a[edge2ind[e]]\n",
    "    return a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Smuggler's objective\n",
    "\n",
    "The smuggler chooses a path $\\mathcal{P}$ only after $c$ is fixed.\n",
    "His goal is to find a valid path in the graph which maximizes $\\sum_{j \\in \\mathcal{P}} c_j$.\n",
    "Note that this problem can be cast and solved exactly as a [shortest path problem in graph theory](http://en.wikipedia.org/wiki/Shortest_path_problem).\n",
    "We formulate it here as a convex problem so that the security team can use the model to thwart the smuggler.\n",
    "\n",
    "We'll denote possible paths with a Boolean decision variable $x \\in \\lbrace 0,1 \\rbrace^m$, where $x_j = 1$ denotes $j \\in \\mathcal{P}$. Note that this allows us to write the smuggler's objective as\n",
    "\n",
    "$$\n",
    "\\sum_{j \\in \\mathcal{P}} c_j = \\sum_{j=0}^{n-1} c_j x_j = c^T x.\n",
    "$$\n",
    "\n",
    "The $x$ variable needs to satisfy certain constraints to represent a valid path. Firstly, the number of times the path exits the source node must be exactly one more than the number of times it enters the source node. Similarly, the number of times the path enters the destination node must be exactly one more than the number of times it exits destination node. For all other nodes in the graph, the number of times the path enters and exits the node must be equal. These constraints can be represented as $Ax = b$, where\n",
    "\n",
    "\\\\[\n",
    "b_i = \\begin{cases} \n",
    "-1 & i = 0 \\\\\n",
    "+1 & i = n-1\\\\\n",
    "0 & \\mbox{otherwise,}\n",
    "\\end{cases}\n",
    "\\\\]\n",
    "\n",
    "\\\\[\n",
    "A_{ij} = \n",
    "\\begin{cases} \n",
    "+1 & \\mbox{if edge $j$ enters node $i$} \\\\\n",
    "-1 & \\mbox{if edge $j$ exits node $i$} \\\\\n",
    "0 & \\mbox{otherwise.}\n",
    "\\end{cases}\n",
    "\\\\]\n",
    "\n",
    "The smuggler's problem can then be written so that $P^\\mathrm{opt}(c)$ is the optimum value of the problem\n",
    "\n",
    "\\\\[\n",
    "\\begin{array}{ll}\n",
    "  \\mbox{maximize} & c^T x \\\\\n",
    "  \\mbox{subject to} & A x = b\\\\\n",
    "  & x \\in \\lbrace 0,1 \\rbrace^m.\n",
    "\\end{array}\n",
    "\\\\]\n",
    "\n",
    "This is a linear problem with Boolean variables, which is not convex, but we can transform it into a convex problem by relaxing the Boolean constraints to get the linear program\n",
    "\n",
    "\\\\[\n",
    "\\begin{array}{ll}\n",
    "  \\mbox{maximize} & c^T x \\\\\n",
    "  \\mbox{subject to} & A x = b\\\\\n",
    "  & x \\geq 0,\n",
    "\\end{array}\n",
    "\\\\]\n",
    "with optimum value $P^{\\mathrm{opt}}_\\mathrm{R}(c)$.\n",
    "\n",
    "It turns out that this relaxed convex problem solves the original Boolean problem *exactly*.\n",
    "That is, $P^{\\mathrm{opt}}(c) = P^{\\mathrm{opt}}_\\mathrm{R}(c)$, so we will only write $P^{\\mathrm{opt}}$ in the remainder.\n",
    "In addition, the solution to the LP, $x^\\star$, will be Boolean when there is a unique optimal path.\n",
    "In the presence of multiple optimal paths, the $x^\\star$ may have fractional entries, but an optimal path can still be recovered.\n",
    "\n",
    "For the purposes of the security team, we'll only need that the Boolean and LP optimal values are equal.\n",
    "Before continuing on to security's perspective, we'll look at an example of the smuggler choosing an optimal path in a given graph.\n",
    "\n",
    "# Smuggler's path example\n",
    "We'll solve the smuggler's problem on an example network to see an example of an optimal path.\n",
    "\n",
    "We'll create an $N \\times N$ grid of points with small perturbations to make the graph irregular, connecting two nodes with two directed edges (one going each way) if the nodes are within distance $\\eta$ of each other. Pairs of edges $(i,j)$ and $(j,i)$ will have identical evasion probabilities, i.e., the evasion probability is the same in either direction between two connected nodes. The evasion probability will be distributed like $p_j = e^{-a_j}$, where $a_j$ is a uniform random variable over the interval $\\left[0,1\\right]$.\n",
    "\n",
    "The smuggler will find an optimal path from node $0$ in the upper-left corner of the graph plot to node $n-1$ in the lower-right corner of the graph plot.\n",
    "\n",
    "We show the graph with the edge evasion probabilities below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/envs/cvxpy/lib/python3.6/site-packages/networkx/drawing/nx_pylab.py:611: MatplotlibDeprecationWarning: isinstance(..., numbers.Number)\n",
      "  if cb.is_numlike(alpha):\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 161384 characters]\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 360,
       "width": 438
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 10\n",
    "G, pos = formGraph(N,.12,1.2,seed=5)\n",
    "n = G.number_of_nodes()\n",
    "\n",
    "a = get_a(G,seed=2)\n",
    "p = np.exp(-a)    \n",
    "\n",
    "fig, ax = showPaths(G,pos,p,Gnodes=[0],Rnodes=[n-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We form the smuggler's relaxed convex problem and solve it to find his optimal path. We plot the path below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----------------------------------------------------------------\n",
      "           OSQP v0.4.1  -  Operator Splitting QP Solver\n",
      "              (c) Bartolomeo Stellato,  Goran Banjac\n",
      "        University of Oxford  -  Stanford University 2018\n",
      "-----------------------------------------------------------------\n",
      "problem:  variables n = 354, constraints m = 808\n",
      "          nnz(P) + nnz(A) = 1416\n",
      "settings: linear system solver = qdldl,\n",
      "          eps_abs = 1.0e-03, eps_rel = 1.0e-03,\n",
      "          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,\n",
      "          rho = 1.00e-01 (adaptive),\n",
      "          sigma = 1.00e-06, alpha = 1.60, max_iter = 4000\n",
      "          check_termination: on (interval 25),\n",
      "          scaling: on, scaled_termination: off\n",
      "          warm start: on, polish: on\n",
      "\n",
      "iter   objective    pri res    dua res    rho        time\n",
      "   1  -8.9458e+02   8.00e+00   9.97e+01   1.00e-01   2.10e-03s\n",
      "  75   5.0606e+00   1.47e-03   5.64e-05   1.00e-01   8.35e-03s\n",
      "\n",
      "status:               solved\n",
      "solution polish:      unsuccessful\n",
      "number of iterations: 75\n",
      "optimal objective:    5.0606\n",
      "run time:             1.15e-02s\n",
      "optimal rho estimate: 5.09e-01\n",
      "\n"
     ]
    }
   ],
   "source": [
    "A = nx.incidence_matrix(G,oriented=True).toarray()\n",
    "n,m = A.shape\n",
    "\n",
    "b = np.zeros(n)\n",
    "b[0] = -1\n",
    "b[n-1] = 1\n",
    "\n",
    "c = np.log(p)\n",
    "\n",
    "edge2ind = {e: i for i,e in enumerate(G.edges())}\n",
    "\n",
    "B = np.zeros((int(m/2),m))\n",
    "count = 0\n",
    "for i in G:\n",
    "    for j in G[i]:\n",
    "        if i < j:\n",
    "            B[count,edge2ind[(i,j)]] = 1\n",
    "            B[count,edge2ind[(j,i)]] = -1\n",
    "            count += 1\n",
    "\n",
    "\n",
    "x = cp.Variable(shape=m)\n",
    "constr = [A*x == b,x>=0, x <= 1]\n",
    "cp.Problem(cp.Maximize(x.T*c),constr).solve(verbose=True)\n",
    "x = np.array(x.value).flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The evasion probability of the smuggler's optimal path is 6.341943e-03, or 0.634%.\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 150944 characters]\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 360,
       "width": 438
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "path = list(np.flatnonzero(x > .1))\n",
    "showPaths(G,pos,p,path,Gnodes=[0],Rnodes=[n-1])\n",
    "print(\"The evasion probability of the smuggler's \"\n",
    "      \"optimal path is %e, or %.3f%%.\"%(np.exp(x.dot(c)), np.exp(x.dot(c))*100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We run a discrete graph-theoretic shortest path algorithm on the same graph to check that we get the same optimal path.\n",
    "The function `optPath` is using the NetworkX Python package and [Dijkstra's_algorithm](http://en.wikipedia.org/wiki/Dijkstra%27s_algorithm) to compute the optimal path."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The evasion probability of the smuggler's optimal path is 6.343747e-03, or 0.634%.\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 150944 characters]\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 360,
       "width": 438
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "y = optPath(G,p)\n",
    "path = list(np.flatnonzero(y > .1))\n",
    "showPaths(G,pos,p,path,Gnodes=[0],Rnodes=[n-1])\n",
    "print(\"The evasion probability of the smuggler's \"\n",
    "      \"optimal path is %e, or %.3f%%.\"%(np.exp(y.dot(c)), np.exp(y.dot(c))*100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Security's objective\n",
    "\n",
    "The security team's goal is to minimize $P^\\mathrm{opt}(c)$, subject to some constraints (say, limits on budget or personnel), which we'll denote by $c \\in \\mathcal{C}$:\n",
    "\n",
    "\\\\[\n",
    "\\begin{array}{ll}\n",
    "  \\mbox{minimize} & P^{\\mathrm{opt}}(c) \\\\\n",
    "  \\mbox{subject to} & c \\in \\mathcal{C}.\n",
    "\\end{array}\n",
    "\\\\]\n",
    "\n",
    "But note that $P^{\\mathrm{opt}}(c)$ is the optimal value of the problem\n",
    "\n",
    "\\\\[\n",
    "\\begin{array}{ll}\n",
    "  \\mbox{maximize} & c^T x \\\\\n",
    "  \\mbox{subject to} & A x = b\\\\\n",
    "  & x \\geq 0.\n",
    "\\end{array}\n",
    "\\\\]\n",
    "\n",
    "We'd like to combine these two optimization problems into a single problem for the security team to solve, but this is problematic as the variables of one problem, $x$, are multiplied with the variables of the other, $c$, which is not a convex objective in general. To get around this, we'll take the dual (Chapter 5 of the [Convex Optimization book](http://web.stanford.edu/~boyd/cvxbook/)) of the smuggler's problem.\n",
    "\n",
    "Let $D^\\mathrm{opt}(c)$ denote the optimal value of the dual of the smuggler's problem, which is\n",
    "\n",
    "\\\\[\n",
    "\\begin{array}{ll}\n",
    "  \\mbox{minimize} & b^T \\lambda \\\\\n",
    "  \\mbox{subject to} & A^T \\lambda \\geq c,\n",
    "\\end{array}\n",
    "\\\\]\n",
    "with dual variable $\\lambda$.\n",
    "\n",
    "Duality theory guarantees that $D^\\mathrm{opt}(c) = P^{\\mathrm{opt}}(c)$, which allows us to write the security team's problem as\n",
    "\n",
    "\\\\[\n",
    "\\begin{array}{ll}\n",
    "  \\mbox{minimize} & D^{\\mathrm{opt}}(c) \\\\\n",
    "  \\mbox{subject to} & c \\in \\mathcal{C},\n",
    "\\end{array}\n",
    "\\\\]\n",
    "which we can rewrite as the single optimization problem\n",
    "\n",
    "\\\\[\n",
    "\\begin{array}{ll}\n",
    "  \\mbox{minimize} & b^T \\lambda \\\\\n",
    "  \\mbox{subject to} & A^T \\lambda \\geq c\\\\\n",
    "  &c \\in \\mathcal{C},\n",
    "\\end{array}\n",
    "\\\\]\n",
    "where $c$ and $\\lambda$ are the optimization variables.\n",
    "\n",
    "We will denote the optimal value of this problem as $P^\\star$. By solving to find $c^\\star$, the security team will have optimally allocated resources to make detection of the smuggler as likely as possible.\n",
    "\n",
    "# Security example\n",
    "\n",
    "We'll consider security's problem on the same network as the last example with the edge evasion probabilities modeled as $p_j = e^{-a_j r_j}$, where $r_j \\in \\mathbf{R}_+$ denotes the effort (say, yearly budget) allocated to edge $j$. We'll assume $a_j \\in \\mathbf{R}_{++}$ are given and represent the cost involved in securing an edge. As in the last example, $a_j$ is a uniform random variable over the interval $\\left[0,1\\right]$.\n",
    "\n",
    "We'll use the same random seed as the last example, so the last example corresponds to the specific allocation $r_j = 1$ for all $j$ in the current model. We'll use this to compare the detection probability of a naive, uniform effort allocation against the optimal allocation.\n",
    "\n",
    "For this example, we'll impose a maximum budget constraint $\\sum_{j=0}^{m-1} r_j = m$, and a uniform spending limit on each edge, $r_j \\leq R$.\n",
    "\n",
    "We'll also constrain that the evasion probability is equal in both directions. That is, edge $(i,j)$ and edge $(j,i)$ have equal evasion probabilities. We'll enforce that constraint with $Br = 0$, for some matrix $B$.\n",
    "\n",
    "The final model is\n",
    "\\\\[\n",
    "\\begin{array}{ll}\n",
    "  \\mbox{minimize} & b^T \\lambda \\\\\n",
    "  \\mbox{subject to} & A^T \\lambda \\geq c \\\\\n",
    "  & c_j = -a_j r_j \\\\\n",
    "  & Br = 0 \\\\\n",
    "  & 0 \\leq r \\leq R\\\\\n",
    "  &\\mathbf{1}^T r = m\n",
    "\\end{array}\n",
    "\\\\]\n",
    "\n",
    "We solve the model below with $R=5$ and report the evasion probability of the smuggler's optimal path."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----------------------------------------------------------------\n",
      "           OSQP v0.4.1  -  Operator Splitting QP Solver\n",
      "              (c) Bartolomeo Stellato,  Goran Banjac\n",
      "        University of Oxford  -  Stanford University 2018\n",
      "-----------------------------------------------------------------\n",
      "problem:  variables n = 454, constraints m = 1240\n",
      "          nnz(P) + nnz(A) = 2478\n",
      "settings: linear system solver = qdldl,\n",
      "          eps_abs = 1.0e-03, eps_rel = 1.0e-03,\n",
      "          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,\n",
      "          rho = 1.00e-01 (adaptive),\n",
      "          sigma = 1.00e-06, alpha = 1.60, max_iter = 4000\n",
      "          check_termination: on (interval 25),\n",
      "          scaling: on, scaled_termination: off\n",
      "          warm start: on, polish: on\n",
      "\n",
      "iter   objective    pri res    dua res    rho        time\n",
      "   1  -3.3084e+01   3.54e+02   3.54e+04   1.00e-01   2.59e-03s\n",
      " 100  -2.2269e+01   3.28e-01   1.58e-03   5.35e-03   9.06e-03s\n",
      "\n",
      "status:               solved\n",
      "solution polish:      unsuccessful\n",
      "number of iterations: 100\n",
      "optimal objective:    -22.2692\n",
      "run time:             1.05e-02s\n",
      "optimal rho estimate: 4.10e-03\n",
      "\n"
     ]
    }
   ],
   "source": [
    "nu = cp.Variable(shape=n)\n",
    "r = cp.Variable(shape=m)\n",
    "constr = [A.T*nu >= -cp.multiply(a,r), cp.sum(r) == m, r >= 0, B*r == 0, r <= 5]\n",
    "cp.Problem(cp.Minimize(nu.T*b),constr).solve(verbose=True)\n",
    "nu = np.array(nu.value).flatten()\n",
    "r = np.array(r.value).flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The evasion probability of the smuggler's optimal path is 2.131035e-10.\n",
      "The smuggler's chance of evasion is 29759913.87 times smaller than with the uniform resource allocation.\n"
     ]
    }
   ],
   "source": [
    "print(\"The evasion probability of the smuggler's optimal path is %e.\"%(np.exp(nu.dot(b)),))\n",
    "print(\"The smuggler's chance of evasion is %.2f times smaller than with the uniform resource allocation.\"%(np.exp(x.dot(c))/np.exp(nu.dot(b))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we plot the resulting edge evasion probabilities from the optimal allocation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<Figure size 576x432 with 2 Axes>,\n",
       " <matplotlib.axes._subplots.AxesSubplot at 0xd19325278>)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 148084 characters]\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 360,
       "width": 438
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "p = np.exp(-a*r)\n",
    "showPaths(G,pos,p,Gnodes=[0],Rnodes=[n-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now solve the smuggler's problem with these optimal resource allocations, but we won't recover a Boolean solution for $x^{\\star}$ because the optimal path is not unique. Note, however, that the optimal evasion probability is the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----------------------------------------------------------------\n",
      "           OSQP v0.4.1  -  Operator Splitting QP Solver\n",
      "              (c) Bartolomeo Stellato,  Goran Banjac\n",
      "        University of Oxford  -  Stanford University 2018\n",
      "-----------------------------------------------------------------\n",
      "problem:  variables n = 354, constraints m = 808\n",
      "          nnz(P) + nnz(A) = 1416\n",
      "settings: linear system solver = qdldl,\n",
      "          eps_abs = 1.0e-03, eps_rel = 1.0e-03,\n",
      "          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,\n",
      "          rho = 1.00e-01 (adaptive),\n",
      "          sigma = 1.00e-06, alpha = 1.60, max_iter = 4000\n",
      "          check_termination: on (interval 25),\n",
      "          scaling: on, scaled_termination: off\n",
      "          warm start: on, polish: on\n",
      "\n",
      "iter   objective    pri res    dua res    rho        time\n",
      "   1  -8.8977e+02   8.00e+00   4.71e+02   1.00e-01   9.14e-04s\n",
      " 200   2.0447e+01   7.16e-03   6.89e-03   1.00e-01   6.66e-03s\n",
      " 325   2.0607e+01   9.79e-04   2.46e-03   1.00e-01   9.95e-03s\n",
      "\n",
      "status:               solved\n",
      "solution polish:      unsuccessful\n",
      "number of iterations: 325\n",
      "optimal objective:    20.6072\n",
      "run time:             1.08e-02s\n",
      "optimal rho estimate: 1.37e-01\n",
      "\n",
      "The evasion probability of the smuggler's optimal path is 1.123036e-09.\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 126012 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 250,
       "width": 373
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "c = np.log(p)\n",
    "x = cp.Variable(shape=m)\n",
    "constr = [A*x == b,x>=0, x <= 1]\n",
    "cp.Problem(cp.Maximize(x.T*c),constr).solve(verbose=True)\n",
    "x = np.array(x.value).flatten()\n",
    "\n",
    "plt.plot(x)\n",
    "print(\"The evasion probability of the smuggler's optimal path is %e.\"%(np.exp(x.dot(c)),))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use [Dijkstra's_algorithm](http://en.wikipedia.org/wiki/Dijkstra%27s_algorithm) again to recover an optimal path for the smuggler in the case that the path is not unique and plot it below. We again check that the detection probability is what we predicted previously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The evasion probability of the smuggler's optimal path is 9.553365e-10.\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 135812 characters]\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 360,
       "width": 438
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = optPath(G,p)\n",
    "path = list(np.flatnonzero(x > .1))\n",
    "showPaths(G,pos,p,path,Gnodes=[0],Rnodes=[n-1])\n",
    "print(\"The evasion probability of the smuggler's optimal path is %e.\"%(np.exp(x.dot(c)),))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Guard placement example\n",
    "\n",
    "We can now go through a more advanced example.\n",
    "\n",
    "In the previous example, security had direct control over the evasion probability on each edge. In this section, we'll consider an example where security can only control the edge probabilities through the placement of guards, each of whom will effect edges which are close to them.\n",
    "\n",
    "In the last example, the security team put a lot of effort on edges near the source and destination nodes. We'll make the team's job a bit harder by allowing the smuggler's path to start at any node on the left side of the graph, and end at any node on the right. This will force the team to allocate resources more evenly over the graph.\n",
    "\n",
    "To make things even more interesting, we'll add 'buildings' in the region covered by the graph. The smuggler and security will have to work around the buildings. That is, no node is located inside a building, no edge crosses through a building, and guards can't see through building walls, restricting their line-of-sight and effect on nearby edges.\n",
    "\n",
    "## Multiple sources and destinations\n",
    "\n",
    "To allow for multiple source and destination nodes without changing the convex formulation, we'll add a dummy source node and a dummy sink node. The source node will have edges going to each node in the graph\n",
    "on the left hand side. The sink node will have edges going to each node on the right. We'll ensure that these dummy edges always have evasion probability 1. The plot below demonstrates the idea. We won't plot the dummy nodes or edges, but we will highlight the 'new' source and sink nodes with light green and red as before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 120468 characters]\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 350,
       "width": 469
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# show dummy source and destination node\n",
    "N = 10\n",
    "G, pos = formGraph(N,.1,1.2,seed=1,dummyNodes=True)\n",
    "n = G.number_of_nodes()\n",
    "fig, ax = showPaths(G,pos,Gnodes=[0],Rnodes=[n-1])\n",
    "ax.axis([-.6,1.6,-.05,1.05]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Guards\n",
    "\n",
    "The security team will select a subset of guards to place in the graph. Each guard will have a profile, which gives his effect on the edge evasion probabilities. The guard's effect will be a function of his distance to the center of each edge. In principle, guards could be placed anywhere in the region occupied by the graph, but we'll only consider guards placed on node positions to make them easy to visualize.\n",
    "\n",
    "An example of a few guards and the resulting evasion probabilities is shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 161580 characters]\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 360,
       "width": 438
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 10\n",
    "k = 5\n",
    "G, pos = formGraph(N,.1,1.2,seed=3,dummyNodes=True)\n",
    "n = G.number_of_nodes()\n",
    "guardIdxs = list(np.random.choice(list(range(N+1,n-N-1)),k,replace=False))\n",
    "\n",
    "visibleNodes = range(1,n-1)\n",
    "Gnodes = range(1,N+1)\n",
    "Rnodes = range(n-2,n-N-2,-1)\n",
    "\n",
    "edgeVals = getGuardEffects(G,pos,guardIdxs,dummyNodes=[0,n-1])\n",
    "edgeProbs = edgeVals.sum(axis=1)\n",
    "edgeProbs = np.exp(edgeProbs)\n",
    "\n",
    "fig, ax = showPaths(G,pos,edgeProbs,visibleNodes=visibleNodes,Gnodes=Gnodes,Rnodes=Rnodes,guards=guardIdxs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Buildings\n",
    "\n",
    "We'll also add buildings, which will modify the graph topology and also restrict the view of the guards."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 152140 characters]\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 360,
       "width": 438
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 10\n",
    "k = 5\n",
    "buildings = [(.2,.8,.3,.1),\n",
    "             (.35,.1,.3,.1),\n",
    "             (.55,.55,.3,.1),\n",
    "             (.15,.35,.3,.1),\n",
    "             (.65,.25,.3,.1)]\n",
    "\n",
    "G, pos = formGraph(N,.11,1.25,buildings=buildings,seed=1,dummyNodes=True)\n",
    "n = G.number_of_nodes()\n",
    "guardIdxs = list(np.random.choice(list(range(N+1,n-N-1)),k,replace=False))\n",
    "\n",
    "visibleNodes = range(1,n-1)\n",
    "Gnodes = range(1,N+1)\n",
    "Rnodes = range(n-2,n-N-2,-1)\n",
    "\n",
    "edgeVals = getGuardEffects(G,pos,guardIdxs,buildings=buildings,dummyNodes=[0,n-1])\n",
    "edgeProbs = edgeVals.sum(axis=1)\n",
    "edgeProbs = np.exp(edgeProbs)\n",
    "\n",
    "fig, ax = showPaths(G,pos,edgeProbs,visibleNodes=visibleNodes,Gnodes=Gnodes,Rnodes=Rnodes,guards=guardIdxs)\n",
    "\n",
    "for x,y,w,h in buildings:\n",
    "    rect = plt.Rectangle((x,y),w,h,fc='y',alpha=.3)\n",
    "    ax.add_patch(rect)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example to solve\n",
    "\n",
    "We'll solve a larger example with buildings, guards, and multiple source and destination nodes. We'll consider a guard at every interior node of the graph, and try to find a limited number, say 10, to secure the graph as well as possible.\n",
    "\n",
    "# TODO\n",
    "- optimization formulation\n",
    "- this one turns out to be a 'true' integer program, unlike the last\n",
    "- use iterative reweighting to get boolean solution\n",
    "- use relaxation to give bounds on distance to true optimal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 17\n",
    "buildings = [(.2,.8,.3,.1),\n",
    "             (.35,.1,.3,.1),\n",
    "             (.55,.55,.3,.1),\n",
    "             (.15,.35,.3,.1),\n",
    "             (.65,.3,.3,.1)]\n",
    "\n",
    "G, pos = formGraph(N,.11,1.25,buildings=buildings,seed=0,dummyNodes=True)\n",
    "n = G.number_of_nodes()\n",
    "guardIdxs = list(range(N+1,n-N-1))\n",
    "\n",
    "visibleNodes = range(1,n-1)\n",
    "Gnodes = range(1,N+1)\n",
    "Rnodes = range(n-2,n-N-2,-1)\n",
    "\n",
    "edgeVals = getGuardEffects(G,pos,guardIdxs,buildings=buildings,dummyNodes=[0,n-1])\n",
    "edgeProbs = edgeVals.sum(axis=1)\n",
    "edgeProbs = np.exp(edgeProbs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 257796 characters]\n",
      "text/plain": [
       "<Figure size 1152x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 469,
       "width": 575
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = showPaths(G,pos,edgeProbs,visibleNodes=visibleNodes,Gnodes=Gnodes,Rnodes=Rnodes,guards=guardIdxs)\n",
    "\n",
    "fig.set_size_inches((16,8))\n",
    "for x,y,w,h in buildings:\n",
    "    rect = plt.Rectangle((x,y),w,h,fc='y',alpha=.3)\n",
    "    ax.add_patch(rect)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----------------------------------------------------------------\n",
      "           OSQP v0.4.1  -  Operator Splitting QP Solver\n",
      "              (c) Bartolomeo Stellato,  Goran Banjac\n",
      "        University of Oxford  -  Stanford University 2018\n",
      "-----------------------------------------------------------------\n",
      "problem:  variables n = 448, constraints m = 1305\n",
      "          nnz(P) + nnz(A) = 27688\n",
      "settings: linear system solver = qdldl,\n",
      "          eps_abs = 1.0e-03, eps_rel = 1.0e-03,\n",
      "          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,\n",
      "          rho = 1.00e-01 (adaptive),\n",
      "          sigma = 1.00e-06, alpha = 1.60, max_iter = 4000\n",
      "          check_termination: on (interval 25),\n",
      "          scaling: on, scaled_termination: off\n",
      "          warm start: on, polish: on\n",
      "\n",
      "iter   objective    pri res    dua res    rho        time\n",
      "   1  -1.5225e-02   1.00e+01   1.24e+03   1.00e-01   1.50e-02s\n",
      " 200   9.8653e+00   5.33e-03   1.96e-03   1.98e-02   1.10e-01s\n",
      "\n",
      "status:               solved\n",
      "solution polish:      unsuccessful\n",
      "number of iterations: 200\n",
      "optimal objective:    9.8653\n",
      "run time:             1.13e-01s\n",
      "optimal rho estimate: 1.08e-02\n",
      "\n",
      "-----------------------------------------------------------------\n",
      "           OSQP v0.4.1  -  Operator Splitting QP Solver\n",
      "              (c) Bartolomeo Stellato,  Goran Banjac\n",
      "        University of Oxford  -  Stanford University 2018\n",
      "-----------------------------------------------------------------\n",
      "problem:  variables n = 448, constraints m = 1305\n",
      "          nnz(P) + nnz(A) = 27688\n",
      "settings: linear system solver = qdldl,\n",
      "          eps_abs = 1.0e-03, eps_rel = 1.0e-03,\n",
      "          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,\n",
      "          rho = 1.00e-01 (adaptive),\n",
      "          sigma = 1.00e-06, alpha = 1.60, max_iter = 4000\n",
      "          check_termination: on (interval 25),\n",
      "          scaling: on, scaled_termination: off\n",
      "          warm start: on, polish: on\n",
      "\n",
      "iter   objective    pri res    dua res    rho        time\n",
      "   1  -5.7726e+02   1.00e+01   1.23e+04   1.00e-01   1.10e-02s\n",
      " 200   1.8639e+01   1.36e-02   4.08e-01   1.00e-01   5.95e-02s\n",
      " 400   1.9265e+01   2.46e-02   7.13e-02   1.92e-02   1.13e-01s\n",
      " 600   1.8962e+01   1.23e-02   5.74e-02   1.92e-02   1.61e-01s\n",
      " 800   1.8931e+01   7.97e-03   4.31e-02   1.92e-02   2.23e-01s\n",
      "1000   1.8566e+01   6.39e-03   3.56e-02   1.92e-02   2.74e-01s\n",
      "1200   1.8689e+01   8.31e-03   1.31e-02   1.92e-02   3.23e-01s\n",
      "1225   1.8471e+01   7.94e-03   3.93e-03   1.92e-02   3.30e-01s\n",
      "\n",
      "status:               solved\n",
      "solution polish:      unsuccessful\n",
      "number of iterations: 1225\n",
      "optimal objective:    18.4712\n",
      "run time:             3.31e-01s\n",
      "optimal rho estimate: 2.86e-02\n",
      "\n"
     ]
    }
   ],
   "source": [
    "A = nx.incidence_matrix(G,oriented=True).toarray()\n",
    "n,m = A.shape\n",
    "numGuards = len(guardIdxs)\n",
    "\n",
    "b = np.zeros(n)\n",
    "b[0] = -1\n",
    "b[n-1] = 1\n",
    "\n",
    "eps = 1e-1\n",
    "w = np.ones(numGuards)\n",
    "\n",
    "for i in range(2):\n",
    "    nu = cp.Variable(shape=n)\n",
    "    r = cp.Variable(shape=numGuards)\n",
    "    constr = [A.T*nu >= edgeVals*r, cp.sum(r) == 10, r >= 0, r <= 1]\n",
    "    cp.Problem(cp.Minimize(nu.T*b/100 + r.T*w),constr).solve(verbose=True)\n",
    "    nu = np.array(nu.value).flatten()\n",
    "    r = np.array(r.value).flatten()\n",
    "    w = 1/(eps+np.abs(r))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0xd1b50d438>]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 15344 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 250,
       "width": 373
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(r,'o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 252016 characters]\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 360,
       "width": 438
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "c = edgeVals.dot(r)\n",
    "edgeProbs = np.exp(c)\n",
    "\n",
    "guards = [guardIdxs[i] for i in range(len(guardIdxs)) if r[i] > .5]\n",
    "\n",
    "fig, ax = showPaths(G,pos,edgeProbs,visibleNodes=visibleNodes,Gnodes=Gnodes,Rnodes=Rnodes,guards=guards)\n",
    "for x,y,w,h in buildings:\n",
    "    rect = plt.Rectangle((x,y),w,h,fc='y',alpha=.3)\n",
    "    ax.add_patch(rect)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The evasion probability of the smuggler's optimal path is 9.213215e-05.\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 252296 characters]\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 360,
       "width": 438
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = optPath(G,edgeProbs)\n",
    "path_inds = list(np.flatnonzero(x > .1))\n",
    "fig, ax = showPaths(G,pos,edgeProbs,path=path_inds,visibleNodes=visibleNodes,Gnodes=Gnodes,Rnodes=Rnodes,guards=guards)\n",
    "print(\"The evasion probability of the smuggler's optimal path is %e.\"%(np.exp(x.dot(c)),))\n",
    "for x,y,w,h in buildings:\n",
    "    rect = plt.Rectangle((x,y),w,h,fc='y',alpha=.3)\n",
    "    ax.add_patch(rect)"
   ]
  }
 ],
 "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
