Network_Flow_Models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MP305 Network Flow Models - Michael Tuite"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview \n",
    "\n",
    "This file contains a number of Python functions for finding the maximal flow through a network $G$ subject to minimal cost using the Ford Fulkerson Algorithm.\n",
    "\n",
    "The network graph $G$ is stored in a list `G` of two element tuples `(i,j)` describing the directed arcs $(i,j)$ of $G$.\n",
    "\n",
    "It is assumed that node number $1$ is the source and the greatest node `Nsink` is the sink.  \n",
    "Thus `G=[(1,2),(2,3),(1,3)]`  describes a network with 3 nodes where node 1 is the source and node 3 is the sink.\n",
    "\n",
    "The capacity, flow and cost of the arc `(i,j)` of `G` are stored in `c[i][j]`, `phi[i][j]` and `l[i][j]`. \n",
    "Here `c`, `phi` and `l` are Python lists.\n",
    "\n",
    "## Python Functions\n",
    "### The `Initialise(G)` function\n",
    "Having defined the network `G`, initialise `c`, `phi` and `l` values to zero via the `Initialise` function before defining their values in any particular example.  The global variable `Nsink`,  the sink node of `G`, is also found by the `Initialise` function .\n",
    "\n",
    "### The main `Iterate(G)`function \n",
    "This implements the full algorithm to find the maximum flow with minimal cost. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The `Iterate(G)` function is based on a number of other Python functions:\n",
    "\n",
    "### `Flows(G)` \n",
    "This checks for conservation of flow and prints out all of the current flows for G and the total cost of this flow. \n",
    "\n",
    "###  `Links(G)`\n",
    "This finds all arcs `(i,j1)`, ` (i,j2)`,  ... *out* of node `i` of `G`.  The nodes `j1,j2,..` are stored in a global list of sets `Out`.\n",
    "\n",
    "###  `SourceSink(G)`\n",
    "This finds all of the paths from source to sink in any network `G` and results in a global set `SinkPaths` of such paths.\n",
    "\n",
    "###  `IncremNet(G)`\n",
    "This finds the Incremental Network `Gp` associated with the current flow of the network `G`.\n",
    "\n",
    "###  `Newflows(G)`\n",
    "This updates the flows `phi` of `G` according to the best chain found through `Gp`. If the maximal flow is reached, then this is indicated and the maximum flow value is outputed. Otherwise, the output is: the change in flow (`eps`), the cost of the best chain, and the best chain.\n",
    "\n",
    "###  `Iterate(G)`\n",
    "Implements the full algorithm to find the maximum flow with minimal cost. \n",
    "The output is as follows:\n",
    "\n",
    "(1) The incremental network `Gp`.\n",
    "\n",
    "(2) The paths through `Gp` from source to sink. \n",
    "\n",
    "(3) The output of `Newflows(G)`. \n",
    "\n",
    "(4) The output of `Flows(G)` giving the current flows and cost of `G`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Initialise (Gin):\n",
    "    global c,phi,l,cp,lp,Nsink\n",
    "    Nsink=1\n",
    "    for arc in Gin:\n",
    "        i,j=arc\n",
    "        Nsink=max(Nsink,i,j)\n",
    "    # for convenience c[i][j] is capacity of arc [i,j]\n",
    "    c=[[0 for i in range(Nsink+1)] for j in range(Nsink+1)]\n",
    "    phi=[[0 for i in range(Nsink+1)] for j in range(Nsink+1)]\n",
    "    l=[[0 for i in range(Nsink+1)] for j in range(Nsink+1)]\n",
    "    cp=[[0 for i in range(Nsink+1)] for j in range(Nsink+1)]\n",
    "    lp=[[0 for i in range(Nsink+1)] for j in range(Nsink+1)]\n",
    "    print(\"All values of c,phi and l initialised to zero\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Flows (Gin):\n",
    "    global Nsink,l,phi\n",
    "    Flowin=[0 for i in range(Nsink+1)]\n",
    "    Flowout=[0 for i in range(Nsink+1)]\n",
    "    for arc in Gin:\n",
    "        i,j=arc\n",
    "        Flowin[j] = Flowin[j] + phi[i][j]\n",
    "        Flowout[i] = Flowout[i] + phi[i][j]\n",
    "    for k in range(2,Nsink):\n",
    "        if Flowin[k] != Flowout[k]:\n",
    "            print(\"*** ERROR *** Flow not conserved at node\", k)\n",
    "    if Flowout[1] != Flowin[Nsink]:\n",
    "        print(\"*** ERROR *** Flow not conserved at source or sink\")\n",
    "    Totalcost = 0\n",
    "    for arc in G:\n",
    "        i,j=arc\n",
    "        phi_ij = phi[i][j]\n",
    "        Totalcost = Totalcost + l[i][j]*phi_ij\n",
    "        print(arc,\" has flow \",phi_ij)\n",
    "    print(\"Total Cost is \", Totalcost)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Links (Gin):\n",
    "    global Nsink,Out\n",
    "    Out=[set() for k in range(Nsink)] # labelled 0..Nsink-1\n",
    "    for arc in Gin:\n",
    "        i,j = arc\n",
    "        Out[i - 1] = Out[i - 1] | set([j])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SourceSink(Gin):\n",
    "# finds all paths SinkPaths from source 1 to sink Nsink of network G \n",
    "    global Nsink,SinkPaths\n",
    "    Links(Gin)\n",
    "    Paths = set() # current paths from source stored as set of tuples\n",
    "    SinkPaths = set() # paths from source to sink Nsink stored as set of tuples\n",
    "    path = 1 # source node label\n",
    "    for node in Out[0]:# need out edge from node 1\n",
    "        pathn = (path,node)\n",
    "        if node == Nsink:\n",
    "            SinkPaths = SinkPaths | set([pathn])\n",
    "        else:\n",
    "            Paths = Paths | set([pathn])\n",
    "    Npaths = len(Paths)\n",
    "    while (0 < Npaths):\n",
    "        NewPaths = set()\n",
    "        for oldpath in Paths:\n",
    "            nold = len(oldpath)\n",
    "            m = oldpath[-1] # last node in tuple oldpath\n",
    "            for mout in Out[m-1]:\n",
    "                if not mout in oldpath:\n",
    "                    if mout == Nsink:\n",
    "                        SinkPaths = SinkPaths | set([oldpath+tuple([Nsink])])\n",
    "                    else:\n",
    "                        NewPaths = NewPaths | set([oldpath+tuple([mout])])\n",
    "        Paths = NewPaths\n",
    "        Npaths = len(Paths)\n",
    "    print(\"Paths from source to sink: \",SinkPaths)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def IncremNet(Gin):\n",
    "# procedure to create incremental network Gp from given flow network G \n",
    "    global Gp,Nsink,phi,c,l,cp,lp,ArcSign\n",
    "# define lists for ArcSign, cp and lp  (indexed by 0..Nsink-1)\n",
    "    cp=[[0 for i in range(Nsink+1)] for j in range(Nsink+1)]\n",
    "    lp=[[0 for i in range(Nsink+1)] for j in range(Nsink+1)]\n",
    "    ArcSign=[[0 for i in range(Nsink+1)] for j in range(Nsink+1)]\n",
    "    Gp=[]\n",
    "    for arc in Gin:\n",
    "        i,j=arc\n",
    "        pij = phi[i][j]; pji = phi[j][i]; cij = c[i][j]; lij = l[i][j]\n",
    "        if (pij < cij and (pji == 0 or not (j,i) in Gin)): # ij arc\n",
    "            #Gp edges, capacitites and costs added\n",
    "            cpij = cij - pij; cp[i][j] = cpij; lpij = lij; lp[i][j] =lpij\n",
    "            ArcSign[i][j] = 1\n",
    "            Gp=Gp + [(i,j)]\n",
    "        if pij>0: # ji arc\n",
    "            cpji = pij; cp[j][i] = cpji; lpji=-lij; lp[j][i] = lpji\n",
    "            ArcSign[j][i] = -1\n",
    "            Gp=Gp + [(j,i)]\n",
    "    print(\"Incremental Network:\",Gp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Newflows(Gin):\n",
    "# A procedure to modify original flows on Gin along SinkPaths of Gp with minimal cost\n",
    "    global Gp,phi,c,l,cp,lp,ArcSign,Out\n",
    "    SourceSink(Gp)\n",
    "    if SinkPaths == set():\n",
    "        Links(Gin)\n",
    "        Flow = 0\n",
    "        for node in Out[0]:\n",
    "            Flow = Flow + phi[1][node]\n",
    "        Cost=0\n",
    "        for arc in Gin:\n",
    "            i,j=arc\n",
    "            Cost=Cost+l[i][j]*phi[i][j]\n",
    "        print(\"Maximal flow found:\", Flow, \" with minimal cost \", Cost)\n",
    "    else:\n",
    "        for k in range(len(SinkPaths)):\n",
    "            cost = 0\n",
    "            epset = set()\n",
    "            path=list(SinkPaths)[k]\n",
    "            for n in range(0, len(path)-1):\n",
    "                i = path[n]; j = path[n+1];  epset = epset | set([cp[i][j]]); cost = lp[i][j] + cost\n",
    "            eps = min(tuple(epset))\n",
    "            if k == 0: # first path\n",
    "                mincost = cost; bestpath = path; besteps = eps\n",
    "            elif cost < mincost: \n",
    "                mincost = cost; bestpath = path; besteps = eps\n",
    "        print(\"A best path in Gp is \", bestpath, \" of minimum cost \", mincost)\n",
    "        print(\"The min capacity on this path is epsilon \", eps)\n",
    "        print(\"The min cost is \", mincost)\n",
    "        for k in range(0, len(bestpath) - 1):\n",
    "            i = bestpath[k]; j = bestpath[k+1]\n",
    "            if ArcSign[i][j] == 1:\n",
    "                phinewij = phi[i][j] + besteps; phi[i][j]=phinewij\n",
    "            else:\n",
    "                phinewji=phi[j][i] = phi[j][i] - besteps; phi[j][i]=phinewji\n",
    "    #print(\"Flow=\",Flow)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Iterate(Gin):\n",
    "    IncremNet(Gin)\n",
    "    Newflows(Gin)\n",
    "    for arc in Gin:\n",
    "        i,j=arc\n",
    "        print((i,j),\" flow = \", phi[i][j])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 1. This is the first example discussed in lectures.\n",
    "\n",
    "## Define `G` and initialise all capacity, flow, cost values to zero  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "G=[(1,2),(1,3),(2,3),(2,4),(2,5),(3,5),(4,6),(5,4),(5,6)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Initialise (G)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Flows(G)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define Capacities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c[1][2] = 4; c[1][3] = 5; c[2][3] = 2  \n",
    "c[2][4] = 1; c[2][5] = 4; c[3][5] = 3  \n",
    "c[4][6] = 2; c[5][4] = 2; c[5][6] = 6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define Flows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phi[1][2]=4; phi[1][3]=1; phi[2][3]=2; phi[2][4]=1; phi[2][5]=1\n",
    "phi[3][5]=3; phi[4][6]=2; phi[5][4]=1; phi[5][6]=3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Flows(G)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use `Iterate` to run one iteration of the Ford Fulkerson Algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Iterate(G)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Iterate(G)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example 2. Power Station/Coal Field example discussed in class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "G2=[(1,2),(1,3),(2,4),(2,5),(3,4),(3,5),(4,6),(5,6)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Initialise(G2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c[1][2]=3; c[1][3]=5; c[2][4]=3; c[2][5]=3 \n",
    "c[3][5]=5; c[3][4]=5; c[4][6]=2; c[5][6]=5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l[1][2]=5; l[1][3]=3; l[2][4]=3 \n",
    "l[2][5]=6; l[3][4]=5 \n",
    "l[3][5]=9 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Flows(G2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Iterate(G2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Continue to iterate until max flow for minimal cost found"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Iterate(G2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Iterate(G2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Iterate(G2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Iterate(G2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}