"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Learning Objectives \n",
"\n",
"* Basic concepts of directed networks\n",
"\n",
"* The min-cost flow problem\n",
"\n",
"* Duality, optimality and equilibrium\n",
"\n",
"### References\n",
"\n",
"* Galichon, *Optimal Transport Methods in Economics*. Ch.8.\n",
"\n",
"* Tolstoi (1930). Methods of finding the minimal total kilometrage in cargo transportation planning in space. *Transportation Planning* [in Russian]\n",
"\n",
"* Koopmans (1949). Optimum utilization of the transportation system. *Econometrica*.\n",
"\n",
"* Schrijver (2002). On the History of the Transportation and Maximum Flow Problems. *Mathematical programming*."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Motivation\n",
"\n",
"### Motivation: network flow problems\n",
"\n",
"In 1930 Tolstoi, a Russian engineer, has to optimize the shipping of cement from factories to cities in the Soviet union through railway.\n",
"\n",
"* each factory produces a fixed number of tons\n",
"\n",
"* each city needs a fixed number of tons -- for now, we'll assume total production=total consumption\n",
"\n",
"* each factory is connected by rail with a few cities, and the corresponding distance is given\n",
"\n",
"* how to ship in order to minimize the total distance travelled?\n",
"\n",
"This problem belongs to the class of *min-cost flow problems*, an important class of linear programming problems, which are the focus of today's lecture. A decade before the invention of linear programming and the work of Kantorovich, Koopmans and Dantzig, Tolstoi described a heuristic method for solving the problem, which led to the optimal solution. As one can recover using modern tools, and will see that his solution was right.\n",
"\n",
"The *shortest path problem*, or how to find the path of minimal distance from a point to another in a network, also belongs in this class. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Fist, let's load libraries we'll need. We'll need `os` and `pandas` for data importation and manipulations, `numpy` and part of `scipy` for computations, `gurobipy` for linear programming using Gurobi, `osmnx` to import map data from [OpenStreetMap](https://www.openstreetmap.org/), `networkx` for operations on networks, parts of `matplotlib` for displays, \n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"ename": "ModuleNotFoundError",
"evalue": "No module named 'osmnx'",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mModuleNotFoundError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[1;32mimport\u001b[0m \u001b[0mosmnx\u001b[0m \u001b[1;32mas\u001b[0m \u001b[0mox\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[1;31mModuleNotFoundError\u001b[0m: No module named 'osmnx'"
]
}
],
"source": [
"import osmnx as ox"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"import os\n",
"import pandas as pd\n",
"import numpy as np\n",
"import scipy.sparse as sp\n",
"\n",
"import networkx as nx\n",
"%matplotlib inline\n",
"import matplotlib as mpl\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib.animation import FuncAnimation\n",
"import gurobipy as grb\n",
"\n",
"thepath = 'https://raw.githubusercontent.com/math-econ-code/mec_optim_2021-01/master/data_mec_optim/'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"# A look at our data\n",
"\n",
"We will be looking at 3 types of data, which we each briefly describe.\n",
"\n",
"\n",
"\n",
"## Soviet Planning Problem\n",
"\n",
"Tolstoi's original data, collected by Schrijver (2002). There are 68 Soviet cities and 10 factories, and 155 links connecting a factory to a city. This yields a sparse $68\\times10$ matrix. Two vectors listing the demand of each city and the supply of each factory are also specified. This is stored in a $69\\times11$ matrix, where we have appended the demand/supply vectors to the right and to the botton of the distance matrix.\n",
"\n",
"We read the data and display it."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
X
\n",
"
Arkhangelsk
\n",
"
Yaroslavl
\n",
"
Murom
\n",
"
Balakhonikha
\n",
"
Dzerzhinsk
\n",
"
Kishert
\n",
"
Sverdlovsk
\n",
"
Artemovsk
\n",
"
Iledzhk
\n",
"
Dekonskaya
\n",
"
demand:
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
Agryz
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
709.0
\n",
"
1064.0
\n",
"
693.0
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
2.0
\n",
"
\n",
"
\n",
"
1
\n",
"
Aleksandrov
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
397.0
\n",
"
NaN
\n",
"
NaN
\n",
"
1180.0
\n",
"
NaN
\n",
"
NaN
\n",
"
4.0
\n",
"
\n",
"
\n",
"
2
\n",
"
Almaznaya
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
81.0
\n",
"
NaN
\n",
"
65.0
\n",
"
1.5
\n",
"
\n",
"
\n",
"
3
\n",
"
Alchevskaya
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
106.0
\n",
"
NaN
\n",
"
114.0
\n",
"
4.0
\n",
"
\n",
"
\n",
"
4
\n",
"
Baku
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
1554.0
\n",
"
NaN
\n",
"
1563.0
\n",
"
10.0
\n",
"
\n",
"
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
...
\n",
"
\n",
"
\n",
"
64
\n",
"
Shchigry
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
566.0
\n",
"
NaN
\n",
"
549.0
\n",
"
4.0
\n",
"
\n",
"
\n",
"
65
\n",
"
Yudino
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
403.0
\n",
"
757.0
\n",
"
999.0
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
0.5
\n",
"
\n",
"
\n",
"
66
\n",
"
Yama
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
44.0
\n",
"
NaN
\n",
"
52.0
\n",
"
5.0
\n",
"
\n",
"
\n",
"
67
\n",
"
Yasinovataya
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
NaN
\n",
"
85.0
\n",
"
NaN
\n",
"
93.0
\n",
"
6.0
\n",
"
\n",
"
\n",
"
68
\n",
"
supply:
\n",
"
5.0
\n",
"
11.5
\n",
"
8.5
\n",
"
12.0
\n",
"
100.0
\n",
"
12.0
\n",
"
15.0
\n",
"
314.0
\n",
"
10.0
\n",
"
55.0
\n",
"
543.0
\n",
"
\n",
" \n",
"
\n",
"
69 rows × 12 columns
\n",
"
"
],
"text/plain": [
" X Arkhangelsk Yaroslavl Murom Balakhonikha Dzerzhinsk \\\n",
"0 Agryz NaN NaN NaN 709.0 1064.0 \n",
"1 Aleksandrov NaN NaN NaN NaN 397.0 \n",
"2 Almaznaya NaN NaN NaN NaN NaN \n",
"3 Alchevskaya NaN NaN NaN NaN NaN \n",
"4 Baku NaN NaN NaN NaN NaN \n",
".. ... ... ... ... ... ... \n",
"64 Shchigry NaN NaN NaN NaN NaN \n",
"65 Yudino NaN NaN NaN 403.0 757.0 \n",
"66 Yama NaN NaN NaN NaN NaN \n",
"67 Yasinovataya NaN NaN NaN NaN NaN \n",
"68 supply: 5.0 11.5 8.5 12.0 100.0 \n",
"\n",
" Kishert Sverdlovsk Artemovsk Iledzhk Dekonskaya demand: \n",
"0 693.0 NaN NaN NaN NaN 2.0 \n",
"1 NaN NaN 1180.0 NaN NaN 4.0 \n",
"2 NaN NaN 81.0 NaN 65.0 1.5 \n",
"3 NaN NaN 106.0 NaN 114.0 4.0 \n",
"4 NaN NaN 1554.0 NaN 1563.0 10.0 \n",
".. ... ... ... ... ... ... \n",
"64 NaN NaN 566.0 NaN 549.0 4.0 \n",
"65 999.0 NaN NaN NaN NaN 0.5 \n",
"66 NaN NaN 44.0 NaN 52.0 5.0 \n",
"67 NaN NaN 85.0 NaN 93.0 6.0 \n",
"68 12.0 15.0 314.0 10.0 55.0 543.0 \n",
"\n",
"[69 rows x 12 columns]"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"soviet_df = pd.read_csv(thepath+'networks_sovietplanning/distances.csv', sep=',',encoding = \"ISO-8859-1\")\n",
"\n",
"soviet_df"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plants and cities are the *nodes* of our network. \n",
"There are *arcs* directly connecting when there is a railway from the plant to the city. In that case, the entry of the corresponding cell is equal to the length of the arc. When there is none, the corresponding cell shows `NaN`.\n",
"We will index the nodes from to `nbPlants+nbCities`, with the plants being indexed by 0,1,...`nbPlants-1`, and the cities being indexed by `nbPlants`, `nbPlants+1`, ... , `nbPlants+nbCities`.\n",
"Let us collect the arcs and their corresponding distances.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"nbPlants = soviet_df.shape[1]-2\n",
"nbCities = soviet_df.shape[0]-1\n",
"arcsList = []\n",
"costslist = []\n",
"nodesList = [soviet_df.columns.to_numpy()[plant+1]+' Plant' for plant in range(nbPlants)] + [soviet_df.iloc[city,0]+' city' for city in range(nbCities)]\n",
"for plant in range(nbPlants):\n",
" for city in range(nbCities):\n",
" if not pd.isna(soviet_df.iloc[city,plant+1]):\n",
" #arcsList.append((plant,nbPlants+city))\n",
" arcsList.append((nodesList[plant],nodesList[nbPlants+city]))\n",
" costslist.append(soviet_df.iloc[city,plant+1])\n",
"nodesList = [soviet_df.columns.values[plant+1]+' Plant' for plant in range(nbPlants)] + [soviet_df.iloc[city,0]+' city' for city in range(nbCities)]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will create a class `network_problem` to collect and story the relevant information we need: nodes, arcs, and arcs length. We will add methods later."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of nodes=78; number of arcs=155.\n"
]
}
],
"source": [
"class network_problem:\n",
" def __init__(self,nodesList,arcsList,costsList):\n",
" self.nbNodes = len(nodesList)\n",
" self.nbArcs = len(arcsList)\n",
" self.nodesList = nodesList\n",
" self.arcsList = arcsList\n",
" self.costsList = costsList\n",
" self.nodesDict = {node:node_ind for (node,node_ind) in zip(self.nodesList,range(self.nbNodes))}\n",
" self.arcsDict = {arc:arc_ind for (arc,arc_ind) in zip(self.arcsList,range(self.nbArcs))}\n",
" print('Number of nodes='+str(self.nbNodes)+'; number of arcs='+str(self.nbArcs)+'.')\n",
"\n",
"soviet_pb = network_problem(nodesList,arcsList,costslist) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## NYC subway network\n",
"\n",
"NYC subway network data: the network data are available on http://www.mta.info. This data is made of two files. The file `nodes.csv` lists the stations: each station is indexed by the line number; each line has the name of the station, and its spatial coordinates. The file `arcs.csv` lists the links between stations: each link specifies the index of the origin, the index of the destination, and the length of the segment. \n",
"\n",
"We shall store this information in a `network_problem` object called `subway_pb`."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of nodes=501; number of arcs=1290.\n"
]
}
],
"source": [
"arcs = pd.read_csv(thepath+'networks_subway/NYC/arcs.csv', sep=',')#.sort_values(by=['route_id'])\n",
"nodes = pd.read_csv(thepath+ 'networks_subway/NYC/nodes.csv', sep=',')\n",
"stopnames = [i+' '+j for i,j in zip(nodes.stop_name,nodes.route_id)]\n",
"arcsList = [(stopnames[o-1],stopnames[d-1]) for (o,d) in list(zip(arcs.from_stop_nb,arcs.to_stop_nb))]\n",
"subway_pb = network_problem(stopnames,arcsList,list(arcs.dis_line)) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Stations' caracteristics are contained in `nodes` dataframe. The relevant information for us will be the name of the station, the horizontal spatial coordinate, and the vertical coordinates."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
stop_name
\n",
"
flow
\n",
"
stop_lat
\n",
"
stop_lon
\n",
"
stop_id
\n",
"
stop_nb
\n",
"
route_id
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
Van Cortlandt Park - 242 St
\n",
"
1000
\n",
"
40.889248
\n",
"
-73.898583
\n",
"
101
\n",
"
1
\n",
"
(1)
\n",
"
\n",
"
\n",
"
1
\n",
"
238 St
\n",
"
1000
\n",
"
40.884667
\n",
"
-73.900870
\n",
"
103
\n",
"
2
\n",
"
(1)
\n",
"
\n",
"
\n",
"
2
\n",
"
231 St
\n",
"
1000
\n",
"
40.878856
\n",
"
-73.904834
\n",
"
104
\n",
"
3
\n",
"
(1)
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" stop_name flow stop_lat stop_lon stop_id stop_nb \\\n",
"0 Van Cortlandt Park - 242 St 1000 40.889248 -73.898583 101 1 \n",
"1 238 St 1000 40.884667 -73.900870 103 2 \n",
"2 231 St 1000 40.878856 -73.904834 104 3 \n",
"\n",
" route_id \n",
"0 (1) \n",
"1 (1) \n",
"2 (1) "
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nodes.head(3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Routes between stations are contained in `arcs` dataframe--the first column is the origin node, the second one is the destination node, and the third one is the length of the arc in meters. "
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
from_stop_nb
\n",
"
to_stop_nb
\n",
"
dis_line
\n",
"
min_time_elapsed
\n",
"
from_stop_id
\n",
"
to_stop_id
\n",
"
route_id
\n",
"
\n",
" \n",
" \n",
"
\n",
"
0
\n",
"
185
\n",
"
10
\n",
"
15000.0
\n",
"
0.18
\n",
"
A09
\n",
"
112
\n",
"
NaN
\n",
"
\n",
"
\n",
"
1
\n",
"
198
\n",
"
23
\n",
"
15000.0
\n",
"
0.18
\n",
"
A24
\n",
"
125
\n",
"
NaN
\n",
"
\n",
"
\n",
"
2
\n",
"
448
\n",
"
25
\n",
"
15000.0
\n",
"
0.18
\n",
"
R16
\n",
"
127
\n",
"
NaN
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" from_stop_nb to_stop_nb dis_line min_time_elapsed from_stop_id \\\n",
"0 185 10 15000.0 0.18 A09 \n",
"1 198 23 15000.0 0.18 A24 \n",
"2 448 25 15000.0 0.18 R16 \n",
"\n",
" to_stop_id route_id \n",
"0 112 NaN \n",
"1 125 NaN \n",
"2 127 NaN "
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"arcs.head(3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Paris streets\n",
"\n",
"Next, we gather data on the network of the Paris streets. To to this, we use the `osmnx` package to retrieve data from [OpenStreetMap](https://www.openstreetmap.org/). Doing this is as simple as specifying an address, and calling the `graph_from_address` function, specifying the desired distance (in meters) around this address. "
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"ename": "NameError",
"evalue": "name 'ox' is not defined",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mNameError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[0mnyu_address\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'57 boulevard Saint-Germain, 75005 Paris, France'\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mstreets_network\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mox\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgraph\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgraph_from_address\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnyu_address\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdist\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m2000\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mnetwork_type\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;34m'walk'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 3\u001b[0m \u001b[0mox\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot_graph\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mox\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mproject_graph\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mstreets_network\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;31mNameError\u001b[0m: name 'ox' is not defined"
]
}
],
"source": [
"nyu_address='57 boulevard Saint-Germain, 75005 Paris, France'\n",
"streets_networkn = ox.graph.graph_from_address(nyu_address, dist=2000,network_type = 'walk')\n",
"ox.plot_graph(ox.project_graph(streets_networkn))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This network is, however, too big for our trial version of Gurobi. We will therefore Zoom in on a smaller area: the small island at the centre of Paris called Île de la Cité."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"address='Île de la cité, Paris, France'\n",
"streets_network = ox.graph.graph_from_address(address, dist=600,network_type = 'walk')\n",
"ox.plot_graph(ox.project_graph(streets_network))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As we see, the essential information we can extract from the graph is the list of nodes (`nodes_roads` below) and the list of arcs (`arcs_roads` below)."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"ename": "NameError",
"evalue": "name 'streets_network' is not defined",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mNameError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0marcslist\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mlist\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mstreets_network\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0medges\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 2\u001b[0m \u001b[0mcostslist\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[0mstreets_network\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m\"length\"\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mj\u001b[0m \u001b[1;32min\u001b[0m \u001b[0marcslist\u001b[0m \u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 3\u001b[0m \u001b[0mstreets_pb\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnetwork_problem\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlist\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mstreets_network\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnodes\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0marcslist\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mcostslist\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;31mNameError\u001b[0m: name 'streets_network' is not defined"
]
}
],
"source": [
"arcslist = [(i,j) for (i,j,k) in list(streets_network.edges)]\n",
"costslist = [streets_network[i][j][0][\"length\"] for i,j in arcslist ]\n",
"streets_pb = network_problem(list(streets_network.nodes),arcslist,costslist) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now that we have the structure of the data in mind, let's develop the conceptual tools we shall need in our analysis."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Topology of networks\n",
"\n",
"### Directed graph\n",
"\n",
"We start by defining the directed graph on which transportation takes place.\n",
"\n",
"---\n",
"**Definition**\n",
"\n",
"A (directed) *graph* $\\left(\\mathcal{Z},\\mathcal{A}\\right)$ is a set of *nodes* (cities) $\\mathcal{Z}$, along with a set of *arcs* $\\mathcal{A}\\subseteq\\mathcal{Z}^{2}$ which are pairs $\\left(x,y\\right)$ where $x\\neq y\\in\\mathcal{Z}$.\n",
"\n",
"---\n",
"Our definition does not allow for an arc to have the same origin and destination. Note that for a dense network, $\\left\\vert \\mathcal{A} \\right\\vert = \\left\\vert \\mathcal{Z}\\right\\vert\\left( \\left\\vert \\mathcal{Z}\\right\\vert - 1 \\right)$. For a line, $\\left\\vert \\mathcal{A}\\right\\vert =\\left\\vert \\mathcal{Z}\\right\\vert -1$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Stylized example**. Take the following example, which we shall use for exposition purposes:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of nodes=4; number of arcs=5.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"C:\\Users\\antoi\\Anaconda3\\lib\\site-packages\\networkx\\drawing\\nx_pylab.py:579: MatplotlibDeprecationWarning: \n",
"The iterable function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use np.iterable instead.\n",
" if not cb.iterable(width):\n",
"C:\\Users\\antoi\\Anaconda3\\lib\\site-packages\\networkx\\drawing\\nx_pylab.py:676: MatplotlibDeprecationWarning: \n",
"The iterable function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use np.iterable instead.\n",
" if cb.iterable(node_size): # many node sizes\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADnCAYAAAC9roUQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxMd/v/8ddkEqIi9rW6SogtfFU1VJuVkNjFTlrctbWC0t5d7tRee6m9i6rWVnvtIWJfG1uKrFUlbag9IrLP74/zQ5UQzMxnluv5eHggM3POe07OXHPNZ875HJ3BYDAghBDCLBxUBxBCCHsiRVcIIcxIiq4QQpiRFF0hhDAjKbpCCGFGUnSFEMKMpOgKIYQZSdEVQggzkqIrhBBmJEVXCCHMSIquEEKYkRRdIYQwIym6QghhRlJ0hRDCjKToCiGEGTmqDiDM71JaJisOJxN3PpXUjBxcnR3xqOBKh1cqU9qlsOp4Qtg0nUxibj+On7vGrB1J7Ey4CEBmTt6d25wdHTAAPtXKMsDbjTrPlVCUUgjbJkXXTiw8cIaxG+PIyMnlYb9xnQ6cHfV8GuRBd68XzZZPCHshwwt2QCu4sdzKznvkfQ0GuJWdy9iNsQBSeIUwMvkizcYdP3eNsRvjClRw/+lWdh5jN8YRk3zNRMmEsE9SdG3crB1JZOTkPtFjM3Jymb0jyciJhLBvUnRt2KW0THYmXHzoGO7DGAywPf4il9MyjRtMCDsmRdeGrTicnO9tf4xvQfbVv+78/9L6qVzd9eN999MBK47kvxwhxOORomvD4s6n3nNY2JPIyMkjLuWGkRIJIaTo2rDUjBwjLSfbKMsRQkjRtWmuzsY5ItDV2ckoyxFCSNG1aR4VXCns+OBfsc6pMIbsu1+Q5d68+sD7OTs64FGxmEnyCWGPpOjasJBXKud7W6FyL3Pz1E4MebncOn2YzHMnHng/AxBSL//lCCEejxRdG1bGpTDeVcui091/W8mAPtxKOsS5aZ25eXIHRdy97ruPTge+1crKJDhCGJHMvWDjjp+7RudvDnAr+/FPkCji5MBPfRriWVkmvxHCWKTTtXF1nivBp0EeFHF6vF+1Ljebyhd/oaaM5wphVFJ07UB3rxf5NKg6RZz0DxxquIchjyJOesJb1ICk3fTo0YOcHOMceiaEkKJrN7p7vUhrl9NUK5pJYUcHnP91VIOzowOFHR3Qp5xkYI1cer1ZlXXr1nHt2jU6d+5MdrYcqyuEMciYrh3IyMhgwIABzJ8/n06dOjHr2wWsOJJMXMoNUjOycXV2wqNiMULqVWZHxHo+//xzoqOj0el0ZGZm0rFjRwCWLVtG4cLypZoQT0OKro2LjY2lZcuWnD17luzsbMLDwxk1alS+98/Ly6N+/fqEh4fTtm1bALKysujatSvp6emsXLmSIkWKmCu+EDZHhhdsWE5ODl5eXpw+ffrO8ICj48PPUnNwcGDUqFF89tln5OVp8zYUKlSIpUuXUqJECVq1akV6errJswthq6To2jBHR0d27dpF9erVcXBwQKfTodfrH/m44OBgihYtyrJly+5Z1o8//sizzz5LUFAQaWlppowuhM2SomvjqlevTnp6OnPmzOHNN9+kSpUqj3yMTqdj9OjRDB8+/J4jF/R6Pd999x3u7u4EBgaSmppqyuhC2CQZ07Vxc+fOZc2aNWzevPmxHmcwGPD29qZ379689dZb99yWl5fHwIEDiY6OZvPmzZQsWdKYkYWwaVJ0bVhGRgbu7u6sXLmSBg0aPPbjd+3axdtvv018fDxOTvfONGYwGBg6dCg7duxg69atlC5d2lixhbBpMrxgw7766ivq1av3RAUXuDMcMX/+/Ptu0+l0TJkyhcDAQHx9ffn777+fNq4QdkE6XRuVnp5OlSpV2Lx5M3Xq1Hni5Rw4cICOHTuSkJCAs7PzfbcbDAZGjBjB8uXL2bZtGxUrVnya2ELYPOl0bdSsWbNo3LjxUxVcAC8vLzw9Pfnmm28eeLtOp2PkyJF069YNb29vkpPlempCPIx0ujboxo0buLm5sX37dmrUqPHUyzty5AgtWrQgKSmJZ555Jt/7TZ48mblz57Jt2zZeeOGFp16vELZIOl0b9OWXX9KkSROjFFyAevXq0bBhQ+bMmfPQ+w0bNoywsDB8fHw4ffq0UdYthK2RTtfGXL16FXd3d/bv34+7u7vRlnvixAn8/f1JSkqiWLGHT/f41VdfMXbsWCIjI6latarRMghhC6TTtTFffPEFrVu3NmrBBahVqxb+/v7MmDHjkfft27cvI0aMwM/Pj9jYWKPmEMLaSadrQy5dukS1atU4fPgwL774otGXHx8fT+PGjUlMTKREiUdfTWLhwoV8+OGHREREULt2baPnEcIaSadrQyZOnEinTp1MUnABqlWrRnBwMFOnTi3Q/bt3787UqVNp0qQJR48eNUkmIayNdLo24vz589SsWZOYmBieffZZk63n9OnTNGjQgPj4+AKfhbZ69Wr69evHunXrnvhEDSFshRRdGzFo0CB0Oh3Tpk0z+br69u1LyZIlGT9+fIEfs379enr16sWaNWto1KiRCdMJYdmk6NqA5ORk6tSpw8mTJ6lQoYLJ13fu3Dnq1KlDbGws5cuXL/DjIiIi6NGjBytWrODNN980YUIhLJcUXRvQv39/XF1dmTBhgtnWGRYWhl6vL/D47m1RUVF07tyZJUuW4O/vb6J0QlguKbpW7vfff6d+/fokJCSYdaavlJQUatasya+//vrYY8i7d++mffv2/PDDDzRr1sxECYWwTFJ0rVyvXr2oXLnyQ697ZioffPAB6enpzJo167Efu3//ftq0acO3335Ly5YtTZBOCMskRdeKJSQk8Prrrxf4uFlju3jxIh4eHhw5cuSJ5lqIjo4mODiY2bNn0759exMkFMLyyHG6VmzkyJEMHjxYScEFKFu2LP369WP06NFP9Pj69esTERHBe++9x9KlS42cTgjLJJ2ulTp58iR+fn4FmgvBlG7P9XDgwAHc3NyeaBknTpygadOmjB8/ntDQUCMnFMKySKdrpYYPH86wYcOUFlyAkiVLEhYWxsiRI594GbVq1SIqKopPPvmEefPmGTGdEJZHOl0rdOzYMYKCgh45v625pKam4ubmxo4dO55qOsnExEQCAgL473//y4ABA4yYUAjLIUXXCrVq1YqAgADCwsJUR7ljwoQJHD58mGXLlj3Vcn7//Xf8/PwYNGgQgwcPNlI6ISyHFF0rc/DgQUJCQkhMTHzgNctUuXnzJlWqVCEiIuKpLxF09uxZ/Pz86NOnDx9++KGREgphGaToWpnAwEDatWtH3759VUe5z7Rp09ixYwdr1qx56mX9+eef+Pv7061bN8LDw42QTgjLIEXXiuzevZvQ0FDi4+MpVKiQ6jj3ycjIwM3NjdWrV/Pqq68+9fLOnz+Pv78/7dq1Y9SoUeh0OiOkFEItKbpWwmAw4Ovry9tvv83bb7+tOk6+5syZw9q1a9m0aZNRlnfx4kUCAgIIDAxkwoQJUniF1ZNDxqxEVFQUKSkpdO/eXXWUh+rduzexsbHs3bvXKMsrW7YsUVFRbNu2jSFDhiA9grB20ulaAYPBQKNGjQgLC6NLly6q4zzSvHnzWLRoEVFRUUZb5rVr12jWrBn16tVj5syZODhIvyCsk+y5VmDTpk3cuHGDTp06qY5SIKGhoZw7d86oRbdEiRJs2bKFX3/9lT59+pCbm2u0ZQthTlJ0LZzBYCA8PJyRI0daTXfn5OTE8OHDCQ8PN+pwgKurK5s2beK3336jZ8+e5OTkGG3ZQpiLdbyK7diaNWswGAy0bdtWdZTH0qVLF65evUpERIRRl+vi4sKGDRs4f/483bt3Jzs726jLF8LUZEzXguXl5VGnTh3GjRtHixYtVMd5bMuXL2fixIkcOnTI6EcdZGRk0L59e5ydnVmyZIlFHkInxINIp2vBli1bRtGiRQkODlYd5Ym0b9+erKws1q5da/RlOzs7s2rVKnJzcwkJCSEzM9Po6xDCFKTTtVA5OTnUqlWLGTNm0KRJE9VxntjatWsJDw/n6NGjJhmTzs7Oplu3bqSmprJ69WqKFCli9HUIYUzS6VqoxYsXU65cOQICAlRHeSotW7akcOHCrFixwiTLd3JyYvHixZQuXZoWLVpw8+ZNk6xHCGORTtcCZWdn4+Hhwfz5823iUuUREREMHjyYEydOoNfrTbKO3Nxc/vOf//Dbb7+xYcMG5fMMC5Ef6XQt0Pfff8/LL79sEwUXoGnTppQuXZrFixebbB16vZ558+ZRvXp1AgMDuX79usnWJcTTkE7XwmRmZuLu7s6yZcvw8vJSHcdotm/fzjvvvENsbCxOTk4mW4/BYCAsLIyDBw8SERFByZIlTbYuIZ6EdLoW5ptvvsHT09OmCi6Ar68vL7zwAj/88INJ16PT6Zg+fTpvvPEGfn5+XLp0yaTrE+JxSadrQdLT03F3d2fdunXUq1dPdRyj27dvH126dCEhIYHChQubdF0Gg4FPP/2UdevWERkZSfny5U26PiEKSjpdCzJnzhy8vLxssuACNGrUiBo1apjl4pM6nY6xY8cSEhKCj48PKSkpJl+nEAUhna6FSEtLw83NjcjISGrVqqU6jslER0fTunVrkpKSzHZM7bhx45g/fz5RUVFUrlzZLOsUIj/S6VqIGTNm4Ovra9MFF6B+/fq8+uqrzJ0712zr/Pjjj+nXrx/e3t6cOXPGbOsV4kGk07UA169fx83NjT179lCtWjXVcUwuJiaGpk2bkpSUhIuLi9nWO3PmTCZPnsy2bduoUqWK2dYrxD9Jp2sBpk6dSnBwsF0UXABPT0+8vb2ZOXOmWdf73nvv8cknn+Dj40N8fLxZ1y3EbdLpKnb58mWqVq3KL7/8wssvv6w6jtnExsbi7e1NYmIixYsXN+u658+fz//+9z+2bt1KjRo1zLpuIaTTVWzy5MmEhITYVcEFqF69Os2aNWPatGlmX3fPnj2ZOHEiAQEBxMTEmH39wr5Jp6vQ33//TfXq1Tl27BjPPfec6jhml5SUhJeXFwkJCZQqVcrs61++fDkDBw5k48aNNnuYnrA80ukqNH78eLp162aXBRfAzc2Ntm3bMmXKFCXr79ChA3PnzqV58+YcPHhQSQZhf6TTVeTPP//E09OTEydOULFiRdVxlPnjjz+oV68ecXFxlC1bVkmGDRs20LNnT1atWkXjxo2VZBD2Q4quIu+++y5FihRh8uTJqqMoZwnbYsuWLXTv3p1ly5bh4+OjLIewfVJ0FbCE7s6S/PXXX9SqVYuTJ08q7fq3b99Op06dWLx4sdVPHi8slxRdBd555x3KlSvH2LFjVUexGO+//z7Z2dnMmDFDaY7du3fTvn17vv/+e4KCgpRmEbZJiq6Z3f7GPjExUeZ6/YfbR3IcPXqU559/XmmWAwcO0KpVK7755htat26tNIuwPXL0gpmNGjWKsLAwKbj/Uq5cOfr06cOYMWNUR8HLy4tNmzbRt29fk13bTdgv6XTN6PZZWElJSbi6uqqOY3Es7ey848eP06xZM6ZMmULXrl1VxxE2QjpdMxoxYgRDhw6VgpuP0qVLM3DgQEaNGqU6CgB16tQhMjKSDz74gAULFqiOI2yEdLpmEhMTQ2BgIElJSRQtWlR1HIt1e8a13bt34+HhoToOAHFxcTRp0oTPPvuMd955R3UcYeWk6JpJmzZt8Pb2ZsiQIaqjWLxx48YRExPDkiVLVEe5IykpCX9/fz788EPeffdd1XGEFZOiawbR0dG0adOGxMREs10twZrdvorG1q1bqV27tuo4d5w5cwY/Pz8GDhwob57iiUnRNYOgoCBatGjBgAEDVEexGlOmTGHv3r2sWrVKdZR7nDt3Dj8/P3r37s1HH32kOo6wQlJ0TWzfvn107dqV+Ph4k18B15bcvjLy2rVreeWVV1THucdff/2Fn58fXbt2JTw8HJ1OpzqSsCJSdE3M39+frl270rt3b9VRrM7MmTPZtGkTGzZsUB3lPhcuXMDf35/WrVszZswYKbyiwOSQMRPavn07f/zxB6GhoaqjWKV33nmHEydOsH//ftVR7lO+fHm2b9/Ohg0b+OCDD5DeRRSUFF0TMRgMpKSkEBMTg5OTk+o4Vqlw4cKcPHnSYq/gW7ZsWaKioti5cyeDBg2SwisKRIYXTCQnJweDwSAF1wiysrLQ6/Xo9XrVUR4oLS2NjIwMSpQogV6vl6EG8VBSdIUQwoxkeEEIIcxIiq4QQpiRo+oANun4cVi8GG7dgs6doVEj1YmsX24u6HTgYCV9Ql6e9re15BVmI3uEsS1fDqGhULgwNGgA3brBjRuqU1mvrCztb73eugqYg4N15RVmI1+kGdvNm/DPWcQ8PWHhQu1v8Xji4uC77yAyEry9oVcvqF1b63odHLTO1xKdOAFr1oCbG7z6KlSpAgbD3e5Xp7v7B+DKFZg3DwoVgpo1oXFjcHZWl1+YlLwVP63c3LsvJrhbcKOitEL7/PN3uzVRcBcuwMSJ2rbduBEqV4bt27Xb9Pq7Bevf2x/g+nXzZv2no0ehXz/IzobNm2HaNO3nOp2W+3bH/s83jAsXtL9Pn4bx42H6dO15CZsknW5B3O5Sbh8n+qBO6+pVSEmBGjW0/+/YAb//DqVLwxdfaC8k6XYLbskSSEyE//wHKlXSfhYbq23LQoXgrbe04Zt/u3wZgoK0N7oRI8Cc1zjLyoLPP4fMTBg3TvtZWhp8/TUsWKB1vG5uWrdeuzZUr64NQ2Vlac8JtO8BXnwR/vhDul0bJZ1uQdzuUkAbPrjdaf39N8yeDc2bQ0AA7Nt39zE+PtCzJ7RqBdWqwYEDSqJbrawsrfBUqqR1jaBtx4kToV49GDVKK7oDB8LevXcfV6iQVpirV7/7+zBXX5GWpu0fgYHa/7OztU8+778Pe/bAyJFQt642/BAeDlOm3M18u1uPj9eKc06OeTILs5OjFwri4kWYNAkOHdLGFlu2hGHDtK7Ezw8mTNC62Ph47aNt8eL3Pj4z8+4L32Cw3LFIS1Kq1N1hAicnrQhPnaoV0h494KeftE8c69fDjz9q279YMe3Ps8+CoyM0aaI93lzbvEQJrTs/f/5u7t9/194gkpKgTh2tw23TBj78UPsUBFqBdXTUHjt7tvam4uIi+4qNkuGFghg7Fi5dgt69teGDuDjtBeXuro3hLVmijTeWKaN1WS+9pL2Ali3TOpxnnoG5c+G551Q/E+tx48bdDvHNN7XCVb48bN0Kq1drb4CVK0NICLRvD7cnhzcYtC/epk+HlSvvfmw3lx9+gHXrtC/DqlSBpk21DPHxcOoUxMRowyTnz8Po0fDGG9rjtm3T9p369WH4cG34SoquTZKi+yg5OdqLfsoUaNhQ+xiYmal92zxrllZ4g4O1F36ZMtpjcnO1L1G2bYOOHcHLS+1zsGZXrmjDB/7+2pvXP+3ZoxXW5s214gbaR/wZM7Sx3cmTtd+XOQ/dyszU3hSio7V9Z+TI+z/53Ha7qP74I0REaIca3n4ewmZJ0X2Ua9e0TrdJk7sviOxs7QSIqlXhn1f2zcvTXkgWOjGLTXjUNj5/HoYMga5dtWEgcxfdh7mdHe6e6LF5s/bFX2CgdtRD1arg4SEdrg2TMd1HcXXVOtjFi+8W3QMH4PXXtRfN7UN7rO3gfWv1z238zwJ886Z2Ykrx4trvpGXL+++v2oOyvPYarFqljf3+9JN2BMyXX8qRLjZMOt2CuH4d3n5bG5s7dUobUpg4UTv8R1iG5GQIC9O+aLtyRRtbbdtWOyRLCAsiRbegbtyAI0fghRe04yiF5Vq9Gs6ehe7d7x4hIISFkKIrhBBmZEEDXkIIYfuk6D6B3Nxc8v59vr9QypZ+Jzk5OWRkZHBDZqezSVJ0H1NGRgbR0dE4WNK34gK9Xs/BgwfJzMxUHeWpOTo6sm3bNtzc3IiOjlYdRxiZVI7HcOnSJSpXrkzlypVVRxEPULFiRZ577jmuq5xlzEiCg4P5+uuvCQoK4oDM22FT5Iu0xxAWFoajoyNffPGF6igiH/3798fV1ZUJEyaojmIUmzZt4q233mLVqlU0btxYdRxhBFJ0C+jcuXPUrVuX2NhYypUrpzqOyEdycjKenp6cOnWKChUqqI5jFJGRkXTp0oVly5bh6+urOo54SlJ0C6hv376ULFmS8ePHq44iHmHw4MEATLs9gbgN2LFjBx07dmThwoU0lfkZrJoU3QI4ffo0DRo0ID4+ntJysL3FO3/+PDVq1CAmJsamxt/37t1L27ZtmT9/PsHBwarjiCckRbcAevbsyfPPP8/IkSNVRxEF9N///pfU1FTmzJmjOopRHTx4kFatWvHVV1/Rpk0b1XHEE5Ci+wjx8fE0btyYxMRESpQooTqOKKBLly5RrVo1oqOjeemll1THMaojR44QFBTEjBkz6NChg+o44jHJIWOPMHLkSIYMGSIF18qUKVOGAQMGMHr0aNVRjK5evXpEREQQFhbGokWLVMcRj0k63Yc4ceIE/v7+/Pbbb7i4uKiOIx7TtWvXcHd3Z9++fbi7u6uOY3QnT56kadOmjBkzhp49e6qOIwpIiu5DtG/fnkaNGjF06FDVUcQTGjNmDHFxcSxcuFB1FJOIj48nICCA8PBw+vTpozqOKAApuvk4cuQILVq0ICkpiWf+fZkYYTVu3LiBm5sbUVFR1KxZU3Uck/jtt9/w9/dn2LBhvPfee6rjiEeQopuPFi1aEBgYyMCBA1VHEU9p0qRJHDp0iOXLl6uOYjJnzpzBz8+Pd999Vz6ZWTgpug9w4MABOnbsSGJiIoXlygNWLz09nSpVqrBp0ybq1q2rOo7JnDt3Dn9/f95++20++eQT1XFEPqToPkCTJk3o0KGDjJHZkOnTpxMZGcnatWtVRzGplJQU/Pz86NSpE8OHD0cnF7i0OFJ0/2XXrl307NmTuLg4nJycVMcRRpKRkYG7uzsrV66kQYMGquOY1IULFwgICKBly5aMHTtWCq+FkaL7DwaDAW9vb3r37s1bb72lOo4wsq+++opVq1YRERGhOorJXbp0iSZNmuDn58fkyZOl8FoQOTniHyIjI7lw4QLdunVTHUWYQM+ePUlISGDPnj2qo5hcmTJl2LZtG7t27SIsLMxmrqphC6TT/f8MBgMNGzZk8ODBdO7cWXUcYSLz589nwYIFbN++3S66v+vXr9O8eXNq1arF3Llz5YonFkB+A//fxo0buXnzJh07dlQdRZhQjx49+Ouvv4iKilIdxSyKFy9OREQE8fHx9O7dm9zcXNWR7J4UXbQuNzw8nFGjRkknYOMcHR0ZMWIE4eHh2MuHvGLFirFx40bOnj1LaGgoOTk5qiPZNakwwOrVq9HpdDJVnp3o1KkTqampbNq0SXUUsylatCjr16/n8uXLdO3alezsbNWR7Jbdj+nm5uZSp04dJkyYIBND25GVK1cybtw4fvnlF7sY270tIyODDh06oNfr+emnn+TkHwXsvtNdtmwZxYoVIygoSHUUYUZt27YlNzeXn3/+WXUUs3J2dmblypXo9XratWtHRkaG6kh2x6473ZycHGrWrMmsWbMICAhQHUeY2fr16/n44485fvy43Y3lZ2dnExoayuXLl1mzZo1M6mRG9rWn/cvChQupWLEi/v7+qqMIBYKDgylatKhNT4STHycnJxYuXEiFChUIDg4mLS1NdSS7YbedblZWFh4eHixYsIA33nhDdRyhyNatWxk4cCAnTpzA0dFRdRyzy83NpW/fvsTFxbFx40ZcXV1VR7J5dtvpzp8/Hzc3Nym4di4gIIBy5cqxePFi1VGU0Ov1fP3113h6etK0aVOuXbumOpLNs8tO9/bkJytWrOC1115THUcotnPnTnr16mXXkxwZDAaGDBnCnj172LJlC6VKlVIdyWbZZaf7zTffULduXSm4AgBvb29efvllvv/+e9VRlNHpdEydOhV/f398fX25ePGi6kg2y+463fT0dNzc3NiwYQP/93//pzqOsBAHDhygU6dOJCQk2PWxqwaDgc8++4xVq1axbds2KlSooDqSzbG7Tnf27Nk0atRICq64h5eXF7Vr1+bbb79VHUUpnU7H6NGj6dy5Mz4+Pvz555+qI9kcu+p07eEiheLJHTlyhJYtW5KUlESRIkVUx1Fu4sSJfP3110RFRfH888+rjmMz7KrTnT59OgEBAVJwxQPVq1eP1157jTlz5qiOYhE+/PBD3nvvPby9vfn9999Vx7EZdtPpXrt2DXd3d/bu3UvVqlVVxxEW6tdff6VJkyYkJSXh4uKiOo5FmDNnDuPHjycyMhJ3d3fVcaye3XS6X3zxBS1atJCCKx6qdu3a+Pr6MmPGDNVRLEb//v357LPP8PX1JTY2VnUcq2cXne6lS5eoVq0a0dHRvPTSS6rjCAsXHx9P48aNSUpKonjx4qrjWIwffviBjz76iC1btlCrVi3VcayWXXS6kyZNokOHDlJwRYFUq1aN4OBgpk6dqjqKRQkNDeWLL76gSZMmHDt2THUcq2Xzne6FCxeoUaMGx48fp3LlyqrjCCtx+vRpGjRoQEJCgpyd9S8rV65kwIABbNiwgfr166uOY3VsvugOGTKEvLw8vvzyS9VRhJXp27cvpUuX5vPPP1cdxeKsW7eO3r178/PPP9OwYUPVcayKTRfd5ORkPD09OXXqlJxZIx7buXPnqFu3LrGxsZQrV051HIuzefNmQkNDWblypUwc9RhsuugOGDAAFxcXJk6cqDqKsFIDBw6kUKFCTJkyRXUUi7Rt2za6dOnC0qVL8fPzUx3HKths0T1z5gyvvPIK8fHxlClTRnUcYaVSUlKoWbMmJ06coFKlSqrjWKSdO3fSoUMHfvzxRwIDA1XHsXg2W3R79+5NpUqVGD16tOoowsoNGzaMjIwMZs6cqTqKxdq3bx9t2rThu+++o0WLFqrjWDSbLLqJiYk0bNiQxMRESpYsqTqOsHIXL17Ew8ODI0eO8MILL6iOY7EOHTpEy5YtmTt3Lm3btlUdx2LZ5HG6I0eOZNCgQVJwhVGULVuWfv36MWbMGNVRLFqDBg3YtGkT/fv356efflIdx2LZXKd76kvwf8kAAAzxSURBVNQpfHx8SEpKkus9CaO5cuUKVatW5eDBg1SpUkV1HIv266+/EhgYyMSJE+nevbvqOBbH5jrd4cOHM2zYMCm4wqhKlSpFWFgYo0aNUh3F4tWuXZvIyEg++ugjvvvuO9VxLI5NdbrHjh2jefPmJCUlUbRoUdVxhI1JTU3Fzc2NXbt24eHhoTqOxUtISCAgIIBPPvmEfv36qY5jMWyq6LZu3Ro/Pz8GDRqkOoqwURMmTODo0aMsXbpUdRSrcPr0afz8/Hj//fcJCwtTHcci2EzRPXToEO3btycxMRFnZ2fVcYSNunnzJm5ubkRERODp6ak6jlX4448/8PPzo3///gwbNkx1HOVspug2a9aMNm3ayMcYYXJTp05l165drF69WnUUq5GcnIy/vz+hoaF8+umnquMoZRNFd8+ePfTo0YP4+HgKFSqkOo6wcbdu3cLd3Z2ff/6ZV155RXUcq5GSkkJAQAAhISGMGDECnU6nOpISNlF0fX196dGjB7169VIdRdiJ2bNns379ejZu3Kg6ilX5+++/CQgIICgoiHHjxtll4bX6Q8aioqJITk4mNDRUdRRhR3r37s2pU6fYv3+/6ihWpVy5cmzfvp0tW7YwdOhQbKDne2xW3ekaDAZef/113n33Xbp166Y6jrAz8+bNY8mSJURGRqqOYnWuXr1Ks2bNePXVV5k+fToODlbf/xWYVT/TzZs3c/36dTp37qw6irBDoaGhnDlzhh07dqiOYnVKlizJ1q1bOXr0KH379iUvL091JLOx2qJrMBgIDw9n5MiR6PV61XGEHXJycmLEiBGEh4fb5cfkp+Xq6kpERASJiYn06tWL3Nxc1ZHMwmqL7s8//0xOTg7t2rVTHUXYsS5dunD58mW2bNmiOopVcnFxYePGjfz555/06NGDnJwc1ZFMzirHdPPy8qhbty5jx46lZcuWquMIO7d8+XImTZrEwYMH7fLbeGO4desW7dq1w8XFhcWLF+Pk5KQ6kslYZae7fPlyihQpIpMlC4vQvn17MjMzWb9+veooVqtIkSKsWbOGrKwsQkJCyMzMVB3JZKyu083NzaVWrVpMmzZNLg0iLMbPP//M8OHDOXLkiF19E29sWVlZdO3alfT0dFauXEmRIkVURzI6q9s7Fi9eTJkyZWjatKnqKELc0apVK5ycnFi1apXqKFatUKFCLF26lBIlStCqVSvS09NVRzI6q+p0s7Oz8fDwYN68efj4+KiOI8Q9Nm/ezNChQ4mJiZEjap5Sbm4uvXr14uzZs6xbtw4XFxfVkYzGqjrdBQsW8NJLL0nBFRYpMDCQkiVLyrSPRqDX65k/fz5ubm40a9aM1NRU1ZGMxmo63czMTKpWrcrSpUtp2LCh6jhCPND27dvp06cPsbGxODo6qo5j9fLy8hg4cCDR0dFs3rzZJq57aDVFd9asWWzYsEEmGBEWz8/Pj44dO5KWloaDgwPvv/++6khWzWAwMHToUHbs2MHWrVspXbq06khPxSqK7q1bt3Bzc2Pt2rUylZ6waJmZmXz00UdMnz4dvV5P/fr12bdvn+pYVs9gMPDxxx+zceNGIiMjKVeunOpIT8wqxnTnzJlDgwYNpOAKi9emTRtmzpxJXl4e2dnZdnNqq6npdDrGjRtH27Zt8fX1JSUlhevXr/P5559b3SnYFj/olJaWxsSJE9m6davqKEI80rRp02jZsiVnz54lMzPTLk5rNRedTsfIkSMpVKgQb775Jk5OTsTFxeHn54eXl5fqeAVm8Z3ujBkz8PHxoXbt2qqjCPFI1apVIyYm5s78zhcuXFCcyPYMHjyYW7duERsbi06nY/78+aojPRaLGNO9lJbJisPJxJ1PJTUjB1dnRzwquBJYtTgNPKuze/duueS1sDqzZs1i8+bNrFu3Lt99vMMrlSntUlh1VKvSqVMnVq5ceWfoplixYly9evWeY6MteXsrLbrHz11j1o4kdiZcBCAz5+6cms6ODmTn5FDy1l/MG9qJOs+VUBVTiCf2qH3cAPhUK8sAbzfZxwvo5MmTfPvttyxdupQrV66QlZXFwoUL6datm1Vsb2VFd+GBM4zdGEdGTi4PS6ADnJ30fBrkQXevF80VT4inVuB9XAfOjrKPPy6DwcDJkycZM2YMISEhZFSubxXbW0nR1XbGWG5lF3y2+CJODnwaVF12SmEVZB83L2va3mb/Iu34uWuM3Rj3WBsH4FZ2HmM3xhGTfM1EyYQwDtnHzeth2zt5di9unTn2wMep2t5mL7qzdiSRkfNkxy5m5OQye0eSkROJB0pNhWefhffeU53E6tjFPn7sGDRsCDVrgqcn/PSTsijWtr3NWnQvpWWyM+FivuMt1w+sIHlmKGe/6MCfX/e97x3KYIDt8Re5nGa7ExxbjPBw8PZWncJi7dmzh6VLl9538kN++/j1/cu5uPrze352ZetXXNn61T0/s5p9/Jln4Icf4ORJ2LwZBg+Ga6brGM+ePcuMGTNIS0u75+cPqymX1k0hN/UiF1eM4uyUEK4fWHHffVRsb7MW3RWHk/O9LftyMjcOr6fCW1N5/v3llO84Csfi5e+7nw5YcST/5Yh/+eEHrROpUwd69CjYYw4fhgsXQOYszteiRYvo3r07L774IosXL75TfPPbx4vW8ObWb4fJy9TmhzXk5ZIet4eiNe5/Y1Oyjz/uflK1Kri7a/+uVAnKlYOLF00W7+DBgwwZMoRKlSrx+eef3ym+D6spZVoORe9alrIhn/H80BUU9wp54P3Mvb3NekZa3PnUew7huIeDA4bcbLIvn0P/THEcS9xfcAEycvL4ce02Tq740oRJrU/RokUZNWrUvdeWOnkSxo6FvXuhTBm4cgUWLYJJk+5fgJsbrFgBeXkwdCj8+CNs2/bAdeXk5PDRRx+RlZVlomdj+fbs2UNubi7Jycm89dZb9O3bl0WLFhF367kH7uOOxctRqEIV0hP241Lbn4w/YtA5Fabws/cff27qffx///vfvXMXPMl+8k+HDkFWFlSpcs+Pc3JyWL9+PVFRUU+dOTExEQcHB27cuMHw4cMZMWIEQ4cO5VbdjvnXlALKyMkjLuXGU2csKLMW3dSM/E+JdCpZiVL+73B9z2IuXfwD55frUdLvPzgWu39GoUIuJXAr5WbKqFbH2dn5/osiRkVBSIj2QgIoVQq6ddP+5Gf2bAgKgueey/cuOp2OKlWqkJ2dbYTk1unw4cN3/q3T6ShfvjwVKlQgNTb/fbxoDW9uxu7CpbY/N0/tfGCXe5sp9/H7Lvr4JPvJbSkpWme8YAH86zJFOp2OMmXK4Ob29M/jxo0bdy6D5ODggIuLCy+//DL7H1JTHkdqhvn2ZbMWXVfnh6+uaE0fitb0IS8zncubZ3Jtx/eUaTn0vvvVqV6VsE4dTRXTdhgM2kGJ//SoDmb/fti9Wyu+aWlaB+PiAuPH37mrXq+nf//+Jg5v2WJjY9m/fz9vvPEGkyZNokGDBgAs/v1ovo95xqMxV6PmkZN6ifSE/VTsMTnf+5p1H3+S/QS0L1uDg2HMGHjA3Ad6vZ7GjRvTuHHjp464fPnyOxcxmDhxIu3atcPBwYGTP+W/vYH7n1c+XJ3Nd/VhsxZdjwquFHY8/8CPA9mXk8m5cRnnyjXQOTqhcyzEg0bHnR0d8KhYzBxxrZ+/P7RtC0OGQOnS2sfGR3Uwixbd/ff330N09D0FV2gGDRpEr169ePXVV+/5+cP2cf0zxSn8fG0ub5yGY/HyOJV58KcJs+/jT7KfZGVpjwkNhQ4dzBDRnw0bNhAYGHjPhT8ftr0B9EVLkHPt/EOXbe7tbdYv0kJeqZzvbYbcbK7t/J5z07uSPKMHeenXKeEdev/9gJB6+S9H/EPNmvDpp9pRCHXqgEymbTQeHh73FVx4+D4O2hBDxpljFK2Z/9CC2ffxJ9lPli2DXbu0N+a6dbU/xx58PKwxlCpViubNm993peVHbe/iXh24vu8nzk7txPWDD75oqLm3t9nPSOvzYzRbYy889DS9/Oh0EFijPHO71zd+MCGMRPZx87K27W32kyPe9XHD2fHJrpTq7KhngI98gSYsm+zj5mVt29vsRbfOcyX4NMiDIk6Pt2rtPGkPPCvLTEzCssk+bl7Wtr31I0aMGGHWNQKelUtQoogT+09fIfcRnwl0OijipJeJQIRVkX3cvKxpeyudTzcm+RqzdySxPf4iOrSDlG+7Pfelb7WyDPBxk3d/YZVkHzcva9jeFnHliMtpmaw4kkxcyg1SM7JxdXbCo2IxQuqpn+VdCGOQfdy8LHl7W0TRFUIIe2HxF6YUQghbIkVXCCHMSIquEEKYkRRdIYQwIym6QghhRlJ0hRDCjKToCiGEGUnRFUIIM5KiK4QQZiRFVwghzEiKrhBCmJEUXSGEMCMpukIIYUZSdIUQwoyk6AohhBn9P+pRLgyH5UN+AAAAAElFTkSuQmCC\n",
"text/plain": [
"