{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Facility location assignment\n", "\n", "Instructor: Frans de Ruiter\n", "\n", "For assignment and data, see https://www.fransderuiter.com/JADS/\n", "\n", "***" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Setup" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "from pyomo.environ import *\n", "import matplotlib.pyplot as plt # optional for plotting\n", "from matplotlib import cm # optional for plotting" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Read data\n", "Here we load the distance table between the cities." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Unnamed: 0AmsterdamAntwerpAthensBarcelonaBerlinBernBrusselsCalaisCologne...PragueRomeRotterdamStrasbourgStuttgartThe HagueTurinVeniceViennaZurich
Distance (km)
AmsterdamNaN016030821639649875209385280...97318358068370356126414491196861
AntwerpNaN16002766146572370446211237...8701660100544659139109012751180687
AthensNaN308227660331225522627302129762562...2198255128262581242830612250199518862449
BarcelonaNaN16391465331201899913141913991539...167914711565107212631589892132719891036
BerlinNaN649723255218990986782936575...354157369780163671211721108666863
\n", "

5 rows × 38 columns

\n", "
" ], "text/plain": [ " Unnamed: 0 Amsterdam Antwerp Athens Barcelona Berlin \\\n", "Distance (km) \n", "Amsterdam NaN 0 160 3082 1639 649 \n", "Antwerp NaN 160 0 2766 1465 723 \n", "Athens NaN 3082 2766 0 3312 2552 \n", "Barcelona NaN 1639 1465 3312 0 1899 \n", "Berlin NaN 649 723 2552 1899 0 \n", "\n", " Bern Brussels Calais Cologne ... Prague Rome Rotterdam \\\n", "Distance (km) ... \n", "Amsterdam 875 209 385 280 ... 973 1835 80 \n", "Antwerp 704 46 211 237 ... 870 1660 100 \n", "Athens 2627 3021 2976 2562 ... 2198 2551 2826 \n", "Barcelona 913 1419 1399 1539 ... 1679 1471 1565 \n", "Berlin 986 782 936 575 ... 354 1573 697 \n", "\n", " Strasbourg Stuttgart The Hague Turin Venice Vienna Zurich \n", "Distance (km) \n", "Amsterdam 683 703 56 1264 1449 1196 861 \n", "Antwerp 544 659 139 1090 1275 1180 687 \n", "Athens 2581 2428 3061 2250 1995 1886 2449 \n", "Barcelona 1072 1263 1589 892 1327 1989 1036 \n", "Berlin 801 636 712 1172 1108 666 863 \n", "\n", "[5 rows x 38 columns]" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data_path = \"https://www.fransderuiter.com/JADS/Facilitylocation/FacilityLocation.xlsx\"\n", "# Create pandas table\n", "distances = pd.read_excel(data_path, sheet_name=0, header=0, skiprows=2, index_col=1)\n", "\n", "# Show table\n", "distances.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Some data processing" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
AmsterdamAntwerpAthensBarcelonaBerlinBernBrusselsCalaisCologneCopenhagen...PragueRomeRotterdamStrasbourgStuttgartThe HagueTurinVeniceViennaZurich
Distance (km)
Amsterdam016030821639649875209385280904...97318358068370356126414491196861
Antwerp16002766146572370446211237861...8701660100544659139109012751180687
Athens3082276603312255226273021297625623414...2198255128262581242830612250199518862449
Barcelona163914653312018999131419139915392230...167914711565107212631589892132719891036
Berlin649723255218990986782936575743...354157369780163671211721108666863
\n", "

5 rows × 37 columns

\n", "
" ], "text/plain": [ " Amsterdam Antwerp Athens Barcelona Berlin Bern Brussels \\\n", "Distance (km) \n", "Amsterdam 0 160 3082 1639 649 875 209 \n", "Antwerp 160 0 2766 1465 723 704 46 \n", "Athens 3082 2766 0 3312 2552 2627 3021 \n", "Barcelona 1639 1465 3312 0 1899 913 1419 \n", "Berlin 649 723 2552 1899 0 986 782 \n", "\n", " Calais Cologne Copenhagen ... Prague Rome Rotterdam \\\n", "Distance (km) ... \n", "Amsterdam 385 280 904 ... 973 1835 80 \n", "Antwerp 211 237 861 ... 870 1660 100 \n", "Athens 2976 2562 3414 ... 2198 2551 2826 \n", "Barcelona 1399 1539 2230 ... 1679 1471 1565 \n", "Berlin 936 575 743 ... 354 1573 697 \n", "\n", " Strasbourg Stuttgart The Hague Turin Venice Vienna Zurich \n", "Distance (km) \n", "Amsterdam 683 703 56 1264 1449 1196 861 \n", "Antwerp 544 659 139 1090 1275 1180 687 \n", "Athens 2581 2428 3061 2250 1995 1886 2449 \n", "Barcelona 1072 1263 1589 892 1327 1989 1036 \n", "Berlin 801 636 712 1172 1108 666 863 \n", "\n", "[5 rows x 37 columns]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#remove 1st column ('unnamed', because there is a comment in first column it was added in the database)\n", "distances = distances.drop(distances.columns[0],1)\n", "\n", "# Show table with removed column\n", "distances.head()" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Index(['Amsterdam', 'Antwerp', 'Athens', 'Barcelona', 'Berlin', 'Bern',\n", " 'Brussels', 'Calais', 'Cologne', 'Copenhagen', 'Edinburgh', 'Frankfurt',\n", " 'Geneva', 'Genoa', 'Hamburg', 'Le Havre', 'Lisbon', 'London',\n", " 'Luxembourg', 'Lyon', 'Madrid', 'Marseille', 'Milan', 'Munich',\n", " 'Naples', 'Nice', 'Paris', 'Prague', 'Rome', 'Rotterdam', 'Strasbourg',\n", " 'Stuttgart', 'The Hague', 'Turin', 'Venice', 'Vienna', 'Zurich'],\n", " dtype='object')\n" ] } ], "source": [ "## Show index names\n", "cities = distances.columns\n", "print(cities)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Model parameters" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
supplydemandfixed cost
Amsterdam4000100150000
Antwerp4000100150000
Athens4000100150000
Barcelona4000100150000
Berlin4000100150000
\n", "
" ], "text/plain": [ " supply demand fixed cost\n", "Amsterdam 4000 100 150000\n", "Antwerp 4000 100 150000\n", "Athens 4000 100 150000\n", "Barcelona 4000 100 150000\n", "Berlin 4000 100 150000" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Fixed cost for opening a DC (same for each city)\n", "fixed_cost = 150000\n", "\n", "# demand equals 100 for each city\n", "demand_per_city = 100\n", "\n", "# Supply capacity of each city is 4000\n", "M = 4000\n", "\n", "# Create pandas dataframe with supply\n", "city_params_dict = {'supply': M*np.ones(cities.size, int), \n", " 'demand': demand_per_city*np.ones(cities.size, int), \n", " 'fixed cost': fixed_cost*np.ones(cities.size, int)}\n", "city_params = pd.DataFrame(data=city_params_dict, index = cities)\n", "city_params.head()" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "150000" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# extract the fixed costs for a city (in this case \"Athens\") as follows:\n", "i_1 = 'Athens'\n", "city_params.loc[i_1,'fixed cost']" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3082" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# extract the fixed costs for a city (in this case \"Athens\" to \"Amsterdam\") as follows:\n", "i_1 = 'Athens'\n", "j_1 = 'Amsterdam'\n", "distances[i_1][j_1]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Model implementation" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "ename": "NameError", "evalue": "name 'XXXX' 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 11\u001b[0m \u001b[1;31m# Objective\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 12\u001b[0m m.value = Objective(\n\u001b[1;32m---> 13\u001b[1;33m \u001b[0mexpr\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mXXXX\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mcities\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 14\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mXXXX\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mcities\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mj\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mcities\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 15\u001b[0m , sense=minimize)\n", "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m(.0)\u001b[0m\n\u001b[0;32m 11\u001b[0m \u001b[1;31m# Objective\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 12\u001b[0m m.value = Objective(\n\u001b[1;32m---> 13\u001b[1;33m \u001b[0mexpr\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mXXXX\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mcities\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 14\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mXXXX\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mcities\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mj\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mcities\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 15\u001b[0m , sense=minimize)\n", "\u001b[1;31mNameError\u001b[0m: name 'XXXX' is not defined" ] } ], "source": [ "# TODO\n", "# write your optimization model in the cell here by replacing the XXXX in the cells\n", "\n", "# Create model\n", "m = ConcreteModel()\n", "\n", "# Variables\n", "m.locations = Var(cities, within=Binary)\n", "m.shipments = Var(cities, cities, within = NonNegativeReals)\n", "\n", "# Objective\n", "m.value = Objective(\n", " expr=sum( XXXX for i in cities)\n", " + sum( XXXX for i in cities for j in cities)\n", " , sense=minimize)\n", " \n", "# Constraints on supplying only when facility is open\n", "m.supply_restriction = ConstraintList()\n", "for i in cities:\n", " m.supply_restriction.add( XXXX )\n", "\n", "# Constraints on demand fulfillment\n", "m.demand_fulfillment = ConstraintList()\n", "for j in cities:\n", " m.demand_fulfillment.add( XXXX )" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# Create model\n", "m = ConcreteModel()\n", "\n", "# Variables\n", "m.locations = Var(cities, within=Binary)\n", "m.shipments = Var(cities, cities, within = NonNegativeReals)\n", "\n", "# Objective\n", "m.value = Objective(\n", " expr=sum(m.shipments[i,j]*distances[i][j] for i in cities for j in cities)\n", " + sum(m.locations[i]*city_params.loc[i,'fixed cost'] for i in cities)\n", " , sense=minimize\n", " )\n", "\n", "# Constraints on demand fulfillment\n", "m.demand_fulfillment = ConstraintList()\n", "for i in cities:\n", " m.demand_fulfillment.add(sum(m.shipments[j,i] for j in cities) == city_params.loc[i,'demand'])\n", " \n", "# Constraints on supply only if DC is open\n", "m.supply_restriction = ConstraintList()\n", "for i in cities:\n", " m.supply_restriction.add(sum(m.shipments[i,j] for j in cities) <= m.locations[i]*city_params.loc[i,'supply'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solve the model" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Status = optimal \n", "\n" ] } ], "source": [ "# Optimize\n", "solver = SolverFactory('cbc') # Take the cbc solver, glpk is very slow\n", "status = solver.solve(m,tee=False,) # Set tee to True to see log of the solver\n", "\n", "# Print the status of the solved mixed integer linear model once it is done\n", "print(\"Status = %s \\n\" % status.solver.termination_condition)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Show the solution\n" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAR4AAAJTCAYAAADe53n5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdabhcVZn28f9NmAkGhYjIC4RRZgKEKAIxCNLSoICiQVGM0ERwQNsGRFsxQDu02GIrDoAiisg8iNgCMoQZMpCQBERoSWhBRBAZAkmA8Lwf1qpkp1JVZ6raVXXO/buuc1XVHldVcp6z9q5976WIwMysTCu1uwFmNvS48JhZ6Vx4zKx0LjxmVjoXHjMrnQuPmZXOhacEkg6UFJKOb/J2h+ftXtPutpj1RVcVnvwL05efif3cz7fz+mOa/BasjVx0QdJlVb8jSyQ9K+lhSVdIOkbSOr3YzgGSLpI0X9JCSS9JekjSeZL27Wn9lZvzdkpzSo1pnwNGAP8NPFs1b1bLW9R9bga2Af7W7oZYW10KPJCfrw1sAuwFHAJ8XdKxEXFx9UqSXg/8Cng38BJwI/AQsATYHDgI+JikUyJicr2dd1XhqfVGcq9mBPDdiJhfcpO6TkS8CDzY7nZY210SEZcVJ0haBTgWOB34laRFEfHrwvyVgauAccDvgIkR8beqbawJfAp4Y8O9R0RX/wDzgQBG9bDctqRK/QTwMvAYcG71esDTeXvVPwuqtnU6cG9efjEwD/gh8KYa+z4wb+P4Xr6nNYDjST22Z4EX8/avAMYVlhuet3sN8CbgPFJPZhEwG/hQb9sCTAcWAGsC3wIeze/rYeAkYOWq5Yv73hD4JfBU3sYtwFvzcq8j9Ub/nLc3G3hPnfctYCJwa37fi4C5wBeAVRrsv8f3DlxW5981gDF9+dx7+Ler7Gd94Iv581sM/B/wn8BaddYbBZxF+v+8OP+/ugIYXWPZb1faDXw8/9u9CMztQ/sObbDMZ/Iy84qfO/CJPH0OsFoP+2k4v6t6PP0laS9ShV4DuJL0n2E70j/aQZLGR8ScvPi3gIOB3YFzgL/k6S8XNvlh4EhgCumXZAmwI3AMcICkMRHx1ACafDHwHmAm6RdqMemXexzwzrzPopHA3cA/gAuBtYAJpL9aL0fE5b3cr4BfA28h/acPUtf7G8Bo4LAa67wRuAv4K6n4rA98ELhB0ltJxX5V0l/KNYEPAVdK2iUiZi/dsaTc9gmkX75LgReAPYFvAuMkvSciXuvne7+E9G/4IeD3wJ2FbVT+jfv6uTdyNrBH3u8C4ADgRODtkt4ZEa8U3vvbgf8hHfL8Lr/39Umf/f6S9o+IKTX28dXcrt8ANwCr96F9jfwY+DKpGL6d9IcEYFJ+/GZELG60gZ7mt73HMtAfeujxkA4nK8scVDXvqDx9Rr2/KHW2uRGwao3pB+f1Tq+a3useD7BBXvYWQFXzBKxbeF35qx+kXsVKhXljgNeAqb1pC+mvZpB6C2tX7WNmnndInX1/p9hWUnc9gGdIv8zFv5r753nnV+3/03n6Lyn8tczv+fQ876hWvPe+fu49/PtVehR/ATao+n/4uzzvXwvTVyf1vhcAY6u2NYrUi5xHocdZ+P/5LLBtH39feuzx5OWuzMudkF+vRfoDG9To1ff1p6u+1eqnfUgnzn4fheNVgIj4KemXahdJu/R2gxHx54h4ucb0q0j/Sf5pYE0GYHHkf/HC9iMi/l5j2X8AJ0WhNxAR00nvbed8bN5bX42IFwrbWQB8Jb88ss6+v1TV1p/nx3WA46Lw1z0ifkf6ZRpdtZ3Pkk5WTorCX8u83X/P8w6vs/9mvXfo2+feyOkR8URhG68CJ+SXxc/x/aRe1ekRMbVqv/OB77Ks51Ht+xHxQI3pzfB4fhyZH9cnfQv+akT8daAbHwqHWpWCclOd+TcDO+efe3uzQUkrkc5FfBTYgfQLNqywyDP9aShARDwh6WbgXZKmk/7y3Eb6672ozmoPRMTCGtP/THr/a5N+QXvjlhrTpuTHnWvMu7+6XRHxkqTngBci4ska6/yF9McAAEnrAVuQ/rOfmI66VvAS6du4ak157/383BtZ4XOMiLmS/g5sJ2nV/Mdr9zx7S0mTa2xnu/y4DSse6k2ldSr/CFH1uimGQuEZkR+fqDO/Mr3HaxcKzgL+hdRF/h/SL1LlP+ck0gnVgXgv8CXSuYr/yNNeknQRqetbXdiqLyOoeDU/Dqszf4Xla2ybiFgg6UWWfZZFzzXYd6N5qxRer5sfNySdt6hnQY1pzXrv0PfPvZFaBRfSubB1SQXx7yx777V6c0XD62yrVd6cHyvnKv9KOnxdWdL6df6g9NpQONSq/Od/U535G1Qt15CkUaSiMw3YKiKOiIiTImJypK/7B3xntYhYEBFfiojNSd3sj5HOwRwJXDDQ7TewsqQ3VE+UNJx0jN+rz6gfKtu9LSLU4GftFu0faPrnvn6d6W8i/R+pHM5W3vs+Pbz3/6rV5D62qVfy1+p75Jf3wNLLMCrXxe0z0H0MhcIzMz+OrzO/Mr14mLUkP9b6a7lFfvxddRdf0pYs+0vRFBHxaET8gvSP/Tiwn6Q1mrmPKu+oMW18fpxZY96A5XMG80nnZGr9ZW+WRv+u1W0a6Oe+wucoaXtSD+f+wjnCu/PjXn3YdqsdSzq38wjpG8uKs/PjFySt2mgDklZrNH8oFJ4bSNdQvFvS/sUZ+eLDXYBZEVEsPJUTiRvX2N78/DhOhZMRkkaw7B+m3yS9uc6J7rVJvY6XWfYL1AqnSFras5C0FnBafvmzFu73DNLhxDnF/RfasZ6knQa4j7r/ri343E+QVOlNVy6+Oz2/LH6OF5MK279JemeNdknSXrkX0lKSVpH0mdzO14DP5pPiFT8lnWfaEbhc0sga21hd0udZdqha06A/xxMRr0o6gvRV5m8kXQH8L+mk3XtIJx4nVq1WORF9hqSxpO7wyxHxrYj43xzKPBCYIekm4A2kb7KeJl0VvNEAmrwZcJukOaSu7eOk80/vyY9fr/WNWpO8RLoI7/78OVWu49mEdKXrFS3aL8D3SX8EPgbsI+n3pD8Y65Euxd8T+B5w3wD2cR+p+Hxc0jDSZxukX6hmf+7TgNmSLmXZdTzbAnfk9wpARCyU9D7SucIbJd1KuqRhMalAjiV9/msDr9A8H8w9MEiFdRTpeqU3kr4cOSYilgsf59+lg0nXSx0IzJd0A/BH0ue4KamH+AZgcsO9D/T7+Hb/0Psrl7cHLiKd9HuZ9B/rPGCzOsv/C+kKzUWseOXy2qS/Cn/K8x8lfe05gnwFcNW2+nIdz3qkTNotpJPWi/PjjVRde0Hh6t0626pcs7FeT21h+SuXTyf90i8mFekv0uDK5Tr7fpo6V9LW+owK894HXJvXf5l08v/u/JlsMZD3nqfvmT/b5ylcudyXz72Hf7/ilctfYtmVy38mXZxa78rlDUjX5/wBWEg6B/RQ/j97GMtfp9TwOrNetq/ys4T0h/Vh0kWjnwBG9GI7B5B6a4/m34FF+f/KL0jnqxqur7wRG+LyV8hbR0Qrz7EMepIuI12bMzIinm53ezrVUDjHY2YdxoXHzErnwmNmpfM5HjMr3aD/Or0V1ltvvRg1alS7m2FW04wZM56OiBWusekkLjz9MGrUKKZPn97uZpjVJOnRdrehJz7HY2alc+Exs9K58JhZ6Vx4zKx0LjxmVjoXHjMrnQuPmZXOhcfMSufCY2alc+Exs9K58JhZ6Vx4zKx0LjxmVjoXHjMrnQuPmZXOhcfMSufCY2alc+Exs9K58JhZ6Vx4zKx0LjxmVjoXHjMrnQuPmZXOhcfMSufCY2al6+jCI+kQSbOqfl6TtH8ftnGqpH0bzJ8o6czmtNjMeqOjhzCOiCuBKyuvJU0CDgeu6836koZFxMktap6Z9VNH93iKJG0FnAx8FBgn6ZrCvDMlTczP50s6WdLtwAcknSfp0DxvN0l3SrpP0lRJa+dNvFnStZIelvStct+Z2dDT0T2eCkmrAL8Cjo+I/5O0WQ+rLIqIPfO6786PqwIXAxMiYpqk1wEL8/KjgZ2BxcAfJX0/Iv5c1YZJwCSAjTfeuEnvzGxo6pYez2nA/RFxUS+Xv7jGtLcAT0TENICIeD4iXs3zboyI5yJiEfAAsEn1yhFxdkSMiYgxI0eO7MdbMLOKju/xSBoPvB/YpTD5VZYvmqtXrfZirU0BUWc3iwvPl9AFn4tZN+voHo+k1wM/A46IiBcKsx4FtpW0mqQRwD692NyDpHM5u+Vtry3JBcasDTr9F+8Y4I3AjyQVp38DuASYDTwMzOxpQxHxsqQJwPclrUE6v1P3a3Yzax1F1Dv6sHrGjBkT06dPb3czzGqSNCMixrS7HY109KGWmQ1OLjxmVjoXHjMrnQuPmZXOhcfMStdxhUfSFEn/VDXtc5LOlXRZu9plZs3TcYUHuBA4rGraYcDPIuLQNrTHzJqsEwvPZcCBklYDkDQKeDPwmKS5edowSadLmiZptqRP5Onjc4/pMkkPSrpA+crDnFo/RdK9kuZI2jpPH5sT6zPz41va8J7NhpSOKzwR8XdgKvDuPOkwUuizeKXjUcBzEbEbsBtwtKRN87ydgc8B2wKbAXsU1ns6InYBfgQcn6c9CIyLiJ1Jt934eq12SZokabqk6U899dQA36XZ0NZxhScrHm4dll8X7QccIWkWcA+wLrBlnjc1Ih6LiNeAWcCownpX5McZhekjgEtzb+oMYLtaDXI63ax5OrXwXAXsI2kXYI2IuLdqvoDPRMTo/LNpRFyf5zVKmi+uMf004OaI2B54Dysm3c2syTqy8ETEAmAKcC4r9nYg3fr02HyDMCRtJWmtfu5uBPB4fj6xn9swsz7oyMKTXQjsBNS6+ddPSDfsujcfIp1F/5P23wK+IekOYFg/t2FmfeB0ej84nW6dzOl0M7MaXHjMrHQuPGZWOhceMytdxxWeBiHRH/ZjWw2HLzaz9ui4wkP9kGit63kaioiTI+KGprTKzJqmEwtPvZDo7ZJOKARDT6nMl/QHSedIul/S9XkUCXoavrhe2NTMWqvjCk+DkOi7SHmssaQhh3eVNC4vsyXwg4jYDniWNADgUoXhiz8bETuRhrVZSOOwKVXbcEjUrEk6rvBktUKi++WfmcC9wNYsC4bOi4hZ+XkxAFpRb/jiRmHT5TgkatY8nTqg31XAd4ohUUmHA9+IiLOKC+ZDsepg6BpV26s3fHElbHpdsxpuZj3ryB5PnZDodcCRkoYDSNpQ0ht7ucl6wxc3M2xqZr3UqT0eSAXnCvIhV0RcL2kb4K58U8EFwEdIPZyGGgxf/BPSYdm9+U6FTwEHN/+tmFmRQ6L94JCodTKHRM3ManDhMbPSdfI5no41Ywak00w2VPiMRHO5x2NmpXPhMbPSdd2hlqR1gRvzyzeRvk6vZBjGRsTLvdjGp4BnI+KC1rTSzBrpusKTs1yjASRNBhZExLd7u76klSPiBy1qnpn1QtcVnnokbQFcFhGVonQSsHJE/Iek24FbgL2AKyStRxpV9Lt53u3AO0lD3Xw8Iu5sz7swGxqG0jme10XEuIj4bo15ioixwAmkYYxXXKCQTl92ZGdm/TGUCk+t8bkqag1tvJxiOh2cTjcbiMFUeF5l+fdTPRTxiw3WrTW0sZm1yGAqPH8lJdBfL2l14IB2N8jMahs0f90jYpGkrwPTgEdIQxybWQdyOr0fnE63TuZ0uplZDS48ZlY6Fx4zK50Lj5mVbsCFR9K6kmbln79Kejw/f1ZSv79ZkjRR0plV06ZI6uiTZmbWswEXnoj4e0SMzhmpHwNn5OejgdcGun0zG3xafag1rM7QwptLulbSDEm3Sdq6rxuW9KOcnbq/Mpxxnv7Pkh6UdLuk70m6Jk+fLOn4wnJz85hcSPpIHtZ4lqSzJA0b6Bs3s/paXXjqDS18NmkgvV2B44Ef1ll/QuEwbhZQPMz693ytwo7AOyTtmK9YPgvYPyL2pBehqjxkzgRgj9xTWwIcXmM5D2Fs1iStvnJ5haGF84B8bwcu1bIbF69WZ/2LI+LTlReSphTmfVDSJNJ72ADYllRIH4mIeXmZC4FJPbRxH2BXYFpuzxrA36oXioizSQWTMWPG+KpLswFodeGpNbTwSqS7/43u70YlbUrqKe0WEf+QdB4pFNroFuz1QqQCfh4RX+xve8ysb0r/Oj0ingfmSfoAgJKd+riZ15HS5s9JWh/YP09/ENiscu6GdAhVMR/YJe9zF2DTPP1G4NDKcMiS3iBpkz62x8z6oF3X8RwOHCXpPuB+4KC+rBwR9wEz87rnAnfk6QuBTwLX5jsLPgk8l1e7HHhDPld0LPBQXucB4MvA9ZJmA78nHbqZWYsMupCopOERsSCPhf4D4OGIOKOZ+3BI1DqZQ6LtcXTu1dxPuofyWW1uj5lVGTT346nIvZum9nDMrLkGY4/HzDqcC4+Zla7UwiPp33PEYXa+Gvmtkj4nac1erHuwpG0LrydKenMT2zZK0oebtT0zq6+0wiNpd+BAYJeI2BHYF/gz8Dmgx8IDHEy6OrliItCUwiNpZdKwNi48ZiUo8+TyBqTROxcDRMTTko4jFY+bJT0dEXtLWhARwwEkHUoqVmcD7yVlsr5MikKMAS6QtBDYHdgb+A7wNHAvsFlEHChpLPBd0lXTC0kjhf5R0kTSSBSrA2uRit82+Ruxnzf7K3gzW6bMwnM9cLKkh4AbSDms70n6PLB3RDxdb8WIuFPS1cA1EXEZgKT9geMjYnohHDouIuZJurCw+oN5+quS9gW+zrKw6u7AjhHxjKTxeXsH1mpDzoVNAth44437/SGYWYmHWhGxgBTGnEQaA/ji3Otohq1ZMRxaMYIUSJ1L+pp9u8K830fEM73ZQXEk0ZEjPZKo2UCUeh1PRCwBpgBTJM0BPlZrscLz6tFA62kUDj0NuDkiDskZrimFeY1GFzWzFinz5PJbJG1ZmDQaeBR4AVi7MP1JSdtIWgk4pDC9erni60bh0BHA4/n5xAZNrN6+mbVImV+nDwd+LumBHMbcFphMOnH8O0k35+VOAq4BbgKeKKx/EXCCpJmSNgfOA36cTwZD/XDot4BvSLoDaHRnwdnAq5Luk/SvA3urZtbIoAmJlhEOrXBI1DqZQ6LlcjjUrEsMmpCow6Fm3WMw9XjMrEu48JhZ6Vx4zKx0TSk8A0md93L7yw3GZ2bdbcAnl6tS54slrQesClwM/BJ4qcY6w/JVzG3VKe0wG2qa0eNZIXUOHMqy1PnNAJIWSDpV0j3A7pJOljQtDyV8dr7+BknHVS4ylHRRYT87SbpJ0sOSjs7LStLpeRtzJE3I08crD12cX59ZyYVJmp/3fTvwAUm75X3dVdlWEz4TM2ugGV+n9zZ1vhYwNyJOBpD0QEScmp+fT+o1/YZ05fKmufe0TmE/OwJvy9uZKem3pHT5aGAnYD3SaKC39qLNi/IQx+RCMykn4L9ZbwWn082aZ8A9nj6kzpeQxraq2FvSPTks+k6WpcZnk+6z8xHS6J8Vv46IhbmQ3QyMBfYELoyIJRHxJHALsFsvmn0xQC5sa0fEnXn6rxq8T6fTzZqkKRcQ9jJ1vqhyPiXfP+eHwJiI+LOkySxLoh8AjCPd+OsrkioFqTrbEdRPpdcbrriikkpvlGo3sxYZcI+nD6nzokoheFrScNI5IXIifaOIuBk4EViHFC4FOEjS6pLWBcYD04BbgQmShkkaSSpYU/P+t5W0mqQRwD61GhER/wBekPS2POmwvr17M+uPZvR4hgPfz4ctrwL/Szrs+hApdf5EROxdXCEinpV0DjCHNKb5tDxrGPDLXCwEnJGXhVRQfgtsDJwWEX+RdCXpPM99pB7QiRHxVwBJl5AO2x4mDXdcz1HAOZJeJPXanmuwrJk1waBJp/dXJdWen58EbBARn220jtPp1sm6IZ0+aEKiA3CApC+SPotHaXyzMDNrgiFfeCLiYvK3XGZWDme1zKx0LjxmVrqmFR5JS3JAdK6k31RddVxr+eWGDJY0WtI/N6Ed5ykNBGhmHaqZPZ6FETE6IrYHngE+1cPyo1h+yODRQJ8Kj9LQw2bWZVp1qHUXsCHUD3IC3wT2yr2kLwCnki4GnCVpgqS1JJ2bg6QzJR2UtzdR0qWSfgNcn7d/Zg6W/hZ4Y6URDYKoUySdIelWSX/IQdErcgD1P1r0mZhZ1vQeg6RhpCuFf5onvY/aQc6TKAwZLOlJUoTi0/n114GbIuLIfNg2VdINeZvFoYffB7wF2AFYH3gAODcvd2adICrAyxExTtJngV+T8mbPAH+SdEZE/L3qfTkkatYkzezxrKE0ysPfgTcAv8/T+xvk3A84KW9zCilmUfmNLw49PK6w/b+QxuOqqBdEBbg6P84B7o+IJ/KtPR4BNqpujEOiZs3T9HM8wCakG4FVzvH0N4gp4P35vNHoiNg4Iv6Q51UPPbzC5deFIOqhEbEDcA7Lh0UX58fXCs8rr33uyKyFmn6OJyKeA44Djpe0CvWDnI2GJAa4DvhM4bzMznV2eStwWN7+BkAlF1YziGpm7deSk8sRMZMU3DwMuJIU1ryPdBhUCXJWDxl8MylRPiufgD4NWAWYnW/WdVqd3V1JCoLOAX5EOpQjIp4l9XLmAFexLIhqZm025EOi/eGQqHWybgiJ+splMyudC4+Zlc6Fx8xK58JjZqVz4TGz0nXVhXKSlpC+Hl8ZmAd8NH9tbmZdpNt6PH1NwJtZB+q2wlPUYwJeaSjjWyRdIukhSd+UdLikqXm5zfNyIyVdnpPs0yTt0cb3ZTboddWhVkUfEvDkaduQekiPAD+JiLE5lf4Z4HPAf5OG0rld0sakuMY2Vft0Ot2sSbqt8FQS8KOAGdRIwANPSqok4J8HpkXEEwCS/kQa6x3SuaJKrmtfUlyjsp/XSVo7Il6oTIiIs4GzIV253Jq3ZzY0dNuhVn8S8NXJ82IqvVJ4VwJ2LyThNywWHTNrrm4rPECfEvC9dT3w6coLSaOb2V4zW15XFh7odQK+t44DxkiaLekB4Jhmt9fMlnE6vR+cTrdO5nS6mVkNLjxmVjoXHjMrnQuPmZWu4y4grAqC/gH4WES81N5WmVkzdWKPpxgEfZmqr7ZzLqsT221mvdTpv8C3AVtIGpWHGv4hcC+wkaQfSZou6X5Jp1RWkPTPkh6UdLuk70m6Jk+fLOn4wnJzJY3Kzz+Sg6OzJJ2Vs2Bm1iIdW3gkrQzsTzrsgjRM8S8iYueIeBT493ytwo7AOyTtmAfxOwvYPyL2BHoc8lPSNsAEYI8cx1gCHF5juUm50E1/6qmnmvEWzYasTiw8lSDodOD/WJZAfzQi7i4s90FJ9wIzSUMTbwtsDTwSEfPyMhf2Yn/7kMZNn5b3uw+wWfVCHsLYrHk67uQyy4KgS+XU+IuF15sCxwO7RcQ/JJ1HGjm0UVj0VZYvtJWRRgX8PCK+OPCmm1lvdGKPpzdeRypEz0lan3RIBvAgsFnl3A3pEKpiPrALgKRdgE3z9BuBQyW9Mc97g6RNWtl4s6GuE3s8PYqI+yTNBO4n3dzrjjx9oaRPAtdKeprlE+qXA0fkw6lpwEN5nQckfRm4Pn9b9grpdhuPlvaGzIaYQRcSlTQ8IhYoHZ/9AHg4Is5o5j4cErVO5pBoexydezX3AyNI33KZWQfpykOtRnLvpqk9HDNrrsHY4zGzDufCY2alc+Exs9J1deGRtCTnq+ZKulTSmn1c/yeStm1V+8ystq4uPPSQZG9E0rCI+JeIeKB1zTOzWrq98BTdBmwBIOkqSTNycn1SZQFJCySdKukeYHdJUySNycPinFcYAvlf2/UmzIaCQfF1eiHJfm2edGREPCNpDVL48/KI+DuwFjA3Ik7O61U2MRrYMPeckLROjX14CGOzJun2Hk+9JPtxku4D7gY2ArbM05eQohPVHiFlvL4v6d2koY+X43S6WfN0e4+nVpJ9PGks9N0j4iVJU1iWRF+Ux1dfTk647wT8Eymn9UHgyFY23Gwo6/bCU8sI4B+56GwNvK2nFSStB7wcEZdL+hNwXovbaDakDcbCcy1wjKTZwB9Jh1s92RD4WeFezr43j1kLdXXhiYjhNaYtZtn9eRouHxHjCy93aWrjzKyubj+5bGZdyIXHzErnwmNmpXPhMbPSDYrCIykk/Vfh9fGSJufnx0g6om2NM7MVDIrCAywG3pevx1lORPw4In7RhjaZWR2DpfC8CpwNrBDuLA5dLGkLSTdIuk/SvZI2z9NPkDRN0uzicMhm1hqDpfBAGlHicEkjGixzAfCDiNgJeDvwhKT9SFmusaSw6K6SxlWv6CGMzZpn0BSeiHge+AVwXK35ktYmJdCvzMsvioiXgP3yz0zgXtIwyFtWr++QqFnzdPWVyzV8l1Q8flZjXr3hjQV8IyI8DI5ZSQZNjwcgIp4BLgGOqjHveeAxSQcDSFot3yr1OuBIScPz9A0rwxmbWWsMqsKT/Rewwrdb2UdJ9+qZDdwJvCkirgd+BdwlaQ5wGbB2KS01G6IG3RDGZfAQxtbJPISxmVkNLjxmVjoXHjMrnQuPmZWuIwtPo9BnP7Y1StLcpjXOzAasIwsPDUKfZtb9OrXwNAp9vkfSPZJm5sDn+nn6ZEnnS7pJ0sOSjq6x7jBJpxcCoZ/I0zeQdGthHPa9Wv0GzYayTi08UD/0eTvwtojYGbgIOLEwb0fgAGB34GRJb65a9yjguYjYDdgNOFrSpsCHgevyGF07AbOqG+OQqFnzdGxWKyKel1QJfS4szPp/wMWSNgBWBeYV5v06IhYCCyXdTEqcF4vIfsCOkg7Nr0eQAqHTgHMlrQJcFRErFJ6IOJvUC2PMmDG+6tJsADq5xwMp9HkUaczziu8DZ0bEDsAnWDZKKEB1Qah+LeAzETE6/2waEddHxK3AOOBx4HzfsdCstTq68NQJfY4gFQiAj1WtcpCk1SWtC4wn9WSKrgOOzT0bJG0laS1JmwB/i4hzSOOve4wtsxbq6MKTVYc+JwOXSroNeLpq2anAb0mjh54WEX+pmv8T4AHg3vwV+1mkw83xwCxJM4H3A//d5PdgZgWDJiSar/NZEBHfbvW+HBK1TuaQqPxUYB8AACAASURBVJlZDR37rVZfRcTkdrfBzHpn0BSeMs2YAap3I1VrqkFyJsCq+FDLzErnwmNmpeuowpNT6ecXXq8s6SlJ1wxgm6dK2rfB/ImSzuzv9s2s7zrtHM+LwPaS1sjRh3ex7GLBfomIk5vSMjNrmo7q8WS/IwU9AT4EXFiZURyOOL+em++3M0rSHySdI+l+SddLWiMvc14lmyVpN0l35iGMp+ZB/gDeLOnanGr/Vjlv02zo6sTCcxFwmKTVSWnze3q53pak4Ym3A54lXYG8lKRVgYuBz+YhjPdlWfh0NDAB2AGYIGmj6o0X0+ngdLrZQHRc4YmI2cAoUm/nf/qw6rxCqnxG3kbRW4AnImJa3s/zEfFqnndjRDwXEYtIkYpNarRr6RDG4CGMzQai087xVFwNfJuUoVq3MP1Vli+WxWT64sLzJcAaVdsUK6bV663bqZ+L2aDQcT2e7Fzg1IiYUzV9Pjk5LmkXYNM+bPNB0rmc3fL6a0tygTFrg478xYuIx6idEL8cOELSLNItLx7qwzZfljQB+H4+8byQdJ7HzEo2aNLpZXI63TqZ0+lmZjW48JhZ6Vx4zKx0LjxmVrquKzw9BUklvVfSSfn5chELM+sMHfl1eg8aBkkj4mrSBYhm1qG6rseTNQqS1rzNhaSj89DF90m6XNKaefp5kr6Xw6OPFAb7M7MW6dbC058g6RURsVsOiP6B5cfq2gDYEzgQ+GatlT2EsVnzdGXh6WeQdHtJt0maAxwObFeYd1VEvBYRDwDr19nn0pDoyJEOiZoNRDee46moFySt5zzg4Ii4T9LEvF5FMSTq27ibtVg3F55zgeciYo6k8b1Yfm3giTx88eEM8M6GZtZ/XVt4GgRJ6/kK6VzQo8AcUiEyszZwSLQfHBK1TuaQqJlZDS48ZlY6Fx4zK50Lj5mVzoXHzEpXWuFpxfDEvdzvTyRtm5/Pl7Refr6glfs1s/rKvI5nwMMTS1q5MBZWr0TEv/RleTNrvbIPtRqlysfmhPjM/PiWPH2ipEsl/Qa4XtIGkm6VNCsPYbxXXm4/SXdJujcvPzxPnyKp4TUNkk7IyfXZkk5pxRs3s2XKLjyNUuUPAuMiYmfgZODrhXm7Ax+LiHcCHwaui4jRwE7ArHz49GVg34jYBZgOfL43DZK0H2n447GkoYx3lTSuxnJOp5s1SamRiYiYLWkUtVPlI4CfS9qSNOLnKoV5v4+IZ/LzacC5OXN1VUTMkvQOYFvgDkkAqwJ39bJZ++Wfmfn1cFIhurWq7WcDZ0O6crmX2zazGtqR1aqXKj8NuDkiDsnFaUph3ouVJxFxa+6RHACcL+l04B+k4vShfrRHwDci4qx+rGtm/dCOr9PrDU88gmUnmyfWW1nSJsDfIuIc4KekIY3vBvaQtEVeZk1JW/WyPdcBRxbOCW0o6Y29fTNm1nel93gapMq/RTrU+jxwU4NNjAdOkPQKsAA4IiKeyvfYuVDSanm5L9OLIY4j4npJ2wB35cO0BcBHgL/17h2ZWV85nd4PTqdbJ3M63cysBhceMytd196BsJ1mzAB1wJ2ZfZRs3co9HjMrXUcWnmYFShvFJYrh0arpNQcENLPm6dRDrQEHShuRNMzhUbP26cgeT9afQOkaki7KYc+LgTUK6yyQdKqke4Ddi70hSR+X9JCkW4A9ynqDZkNVJxee/gRKjwVeiogdga8BuxbWWQuYGxFvjYjbKxMlbQCcQio47yJlvlZQDImCQ6JmA9Gph1r9DZSOA75XWH92YZ0lwOU1dvVWYEpEPAWQe0orxC2KIVHJIVGzgejkHg8sC5ReWDW9EijdHngPsHphXr2isCgiltSZ50JiVqJOLzx9DZTeShqeGEnbkw7RenIPMF7SuvlWGx8YUIvNrEcdXXgi4rGIqBco/YakO4Bhhek/AobnQ6wTgam92McTwGTS/XtuAO4daLvNrDGHRPshneNpf0jU/3RWi0Oig9Suu6Zf+nb/mHUrFx4zK50Lj5mVzoXHzErnwmNmpevqwuNhiM26U1cXHjPrToOq8EhaW9K8fAUykl4nab6kVSSNlnR3Tq5fKen1eZkpkv5T0tScUN+rve/CbPAbVIUnIl4gDQRYuZ3GYcDlEfEK8AvgCzm5Pgf4amHVlSNiLPC5qulLeQhjs+YZVIUn+wnw8fz848DPJI0A1omIW/L0n5OS7BVX5McZwKhaG42IsyNiTESMGTlyZPNbbTaEDLrCExF3AKPyeOrDImJuL1ZbnB+X0MG3CjEbLAZd4cl+QbqVxs8AIuI54B+F8zcfBW6ps66ZtVi3F541JT1W+Pl8nn4B8HqWv4/Px4DTc3J9NHBqyW01s6yrDysiol7h3BO4LCKeLSw7C3hbjW2MLzx/mjrneMysebq68NQi6fvA/sA/t7stZlbboCs8EfGZdrfBzBrr9nM8ZtaFXHjMrHS9KjydEMaUdJ6kQ9vdDjMbuCHf45E0rOelzKyZ+l14qnsglV6RpEMk3aBkgxy8fJOkYZJOlzQtBzU/kZcfL+kWSZfkZb8p6fAc2pwjafPCbveVdFte7sC8/uqSfpaXnSlp7zx9oqQzC+27RtL4SlurhjP+Z0kPSrpd0vckXdPfz8XMetb0b7Ui4kpJ7wc+Bbwb+GpE/FXSJOC5iNhN0mrAHZKuz6vtBGwDPAM8AvwkIsZK+izwGVJ4E9I1Nu8ANgdulrRF3g8RsYOkrYHrJa0wEmiVynDGJ+chkh8mDYk8T1L14IFACokCkwA23njjvn4sZlbQqkOtzwBfBBZHROUXeT/gCEmzSIPorQtsmedNi4gnImIx8CegUpDmsPwFfZdExGsR8TCpQG1NuljwfICIeBB4lBpDEFcpDme8NfBIRMzLr2sWHodEzZpnID2eV8mFS5KAVQvzNgReA9aXtFJEvAYI+ExEXFfcSD78WVyY9Frh9WtVbawe1CXydhu2LysOc1wczrje+mbWIgPp8cwHds3PDwIqN99amRTO/DDwB6CSn7oOOLZwk66tJK3Vx31+QNJK+bzPZsAfWX7Y4q2AjfP0+cDovPxGwNg623wQ2EzSqPx6Qh/bZGZ91Nsez5qSHiu8/g5wDvBrSVOBG4EX87wvAbdFxG35sGqapN+S7pMzCrg395CeAg7uY3v/SEqVrw8cExGLJP0Q+LGkOaRezsSIWJyHN55HOlybS52hiSNioaRPAtdKeppeDHtsZgPjIYwBScMjYkEuiD8AHo6IM+otP2bMmJg+vf1DGJvV4iGMu8fRuXd2PzACOKvN7TEb1AZdSLQ/cu+mbg/HzJrLPR4zK50Lj5mVrqMKTyvCqA6XmnWejio8ZjY0dHzhkbSJpBtzsPRGSRvn6eflQOedkh6p9GpyOPVMSQ/k64feWNjWPjlIOkfSuTkzhtJoo6dIujfP27otb9ZsiOj4wgOcCfwijwB6AfC9wrwNSFmtA4Fv5mmHAG8BdgCOBt4OKcUOnAdMiIgdSN/oHVvY1tMRsQvwI+D46kbII4maNU03FJ7dgV/l5+eTCk3FVTk0+gDpamZII4ReGBFLIuIvwE15+luAeRHxUH7dp9FEHRI1a55uKDzVipdaF8OlqrNMrfm1eDRRs5J0Q+G5EzgsPz8cuL2H5W8FDss3HtsA2DtPf5A0tPEW+bVHEzVrk077y14rjHoccK6kE0jB0o/3sI0rgXeSwqEPkYtLDpR+HLg0J+inAT9ucvvNrBccEu0Hh0StkzkkamZWgwuPmZXOhcfMSufCY2alc+Exs9J1dOGplVaXdIykIxqsM1nSCpEHM+scnXYdT48iwtfemHW5ju7x1FLs0Ug6LqfQZ0u6qLDYTpJukvSwpKPzslIaQnluTqBPyNPHS5oi6bI8jPEF+abvZtYiXdfjqXISsGkezmadwvQdgbeRhiqemW+PsTswmjRc8nqkYXduzcvvDGwH/AW4A9iDqmiGhzA2a56u6/FUmQ1cIOkjpDG1Kn4dEQsj4mngZtJgfnuyLLX+JClKsVtefmpEPJZHPJ2F0+lmLdXthecA0jhYuwIzcgYL+jbUMSyfcnc63azFurbwSFoJ2CgibgZOBNYBhufZB0laXdK6wHhSIPRWYEJOrY8k3YvHo4aatUGn/2WvlVavGAb8UtIIUm/mjIh4Np8Xngr8ljSO+mkR8RdJV5LO89xH6gGdGBF/9W1OzcrndHo/OJ1unczpdDOzGlx4zKx0LjxmVjoXHjMrXamFp79DFOdYwzVV0zw0sVmXGvI9HknD2t0Gs6Gm7YVH0khJl0ualn/26Mc2Ts7rzpV0dg6EbiNpamGZUZJm5+fz8zq3Ax+QtLmkayXNkHSbr+0xa61OuIDwv0kX/92ex0W/DtimxnJ7SZpVeL0xUDn8OjMiTgWQdD5wYET8RtKqkjaLiEeACcAlhfUXRcSeeZ0bgWMi4mFJbwV+SBoiZymHRM2apxMKz77AtoU7UbxO0toR8ULVcrdFxIGVF5LOK8zbW9KJwJrAG4D7gd+QCs0HSeOqT8g/FRfn7Qwnja9+aaENq1U3MiLOBs6GdAFhn9+lmS3VCYVnJWD3iFjYn5UlrU7qoYyJiD9LmgysnmdfTCooVwAREQ8XVn2xsP9nI2J0v1pvZn3W9nM8wPXApysvJPW1AFSKzNO597L0m66I+BMpbf4Vcg+nWkQ8D8yT9IG8f0naqY9tMLM+KLvwrCnpscLP50lDFI/JdxF8ADimLxuMiGeBc0hDFl9FSqIXXQx8hOXP71Q7HDhK0n2kw7SD+tIGM+sbh0T7wSFR62QOiZqZ1eDCY2alc+Exs9K58JhZ6Vx4zKx0pRSe6lS6pImSzmzStudLWq8Z2zKzcgzpHo+T6Wbt0fbCI+k9ku6RNFPSDZLWz9MnS/q5pOtzr+Z9kr6Vhx++VtIqhc2cIGlq/tkir7/c/Xoqva58b5+bJf2KdNEhkr6Shy/+vaQLlYdINrPWKKvwrCFpVuUHOLUw73bgbRGxM3ARaYysis1Jg/YdBPwSuDkidgAW5ukVz0fEWOBM4Lu9aM9Y4N8jYltJY4D3k4Yxfh9Q88IrSZMkTZc0/amnnurFLsysnrJCoguLIUxJE1n2C/7/gIslbQCsCswrrPe7iHhF0hzSOFrX5ulzWH6Y4QsLj2f0oj1TI6Kynz3JQx7ntv2m1gpOp5s1T9sPtYDvk+6nswPwCZaFPiEPLZzHNH8lluU7XmP5ohk1nr9Kfn9K97tYtbDMi4XnjYY2NrMW6ITCMwJ4PD//WD+3MaHweFd+Pp80pjqkQ7VVqO124D15yOPhLH8IZ2Yt0An345lMumfO48DdwKb92MZqku4hFdIP5WnnAL/Otz+9keV7OUtFxDRJV5OGNn4UmA481482mFkvOZ1OugthRCyQtCZwKzApIu6tt7zT6dbJuiGd3gk9nk5wtqRtSeeXft6o6JjZwLnwABHx4Xa3wWwo6YSTy2Y2xLjwmFnpOrrwSFpf0q8kPZIH27tL0iHtbpeZDUzHFp580d9VwK0RsVlE7AocRrrS2cy6WMcWHtJIni9HxI8rEyLi0Yj4vqRhkk7PwxbPlvQJWBoAnSLpshz6vCAXMCTtk4OocySdK2m1PH2F4Y/b8m7NhpBOLjzbAfW+1j4KeC4idgN2A46WVLnwcGfgc8C2wGbAHnnQv/OACTmasTJwbF7+zIjYLSK2B9YADqQGh0TNmqeTC89yJP1A0n2SpgH7AUfkpPs9wLrAlnnRqRHxWM53zSKFSd8CzIuIh/IyPwfG5ed759tyzCH1srartf+IODsixkTEmJEjR7biLZoNGZ18Hc/9pNtVABARn8p3GpwO/B/wmYi4rriCpPHkYGm2hPQeax4+9TD8sZm1SCf3eG4CVpd0bGHamvnxOuDYys3AJG0laa0G23oQGFW5SRjwUeAWGgx/bGat07E9nogISQcDZ0g6EXiKFPT8AnAp6RDq3nwy+Cng4AbbWiTp46Qw6sqkYY5/HBGLJVWGP57PisMfm1kLOCTaDw6JWifrhpBoJx9qmdkg5cJjZqVz4TGz0rnwmFnpXHjMrHQdUXicQjcbWtpeeJxCNxt62l54aH4KfVdJt+Se03WSNpC0TR5tgrzMKEmz83On081K1gmFp5kp9FVIAwQemntO5wJfi4g/AKtK2iyvOwG4JD93Ot2sZB0XmZD0A9Kwwi+TxrnaUVIlQzWClEJ/mZxCz+tUUujPAtsDv88dl2HAE3ndS4APAt8kFZ7KIIB750jGmsAbSOHUFYYx9hDGZs3TCYWn2Sn0+yNi9xr7uZiU1boi7SYedjrdrD064VCrmSn0PwIjJe2el19F0nYAEfEnUoH6CqkIgdPpZm3R9h5Pk1PoL+fDsu9JGkF6f98l9aogFZzTycMkR8SzTqeblc/p9H5wOt06mdPpZmY1uPCYWelceMysdC48Zla6UgqPpCWSZhV+RjVpu1MkrXASTdIHJP1B0s193N7nJK3Z85JmNhBlfZ2+MCJG15spaeWIeLWJ+zsK+GRE9LrwSBpGimD8EnipiW0xsyptu45H0kTgANJFfGtJei/wa+D1wCrAlyPi17l39DvgduDtwOPAQRGxsLCtlYCfAX8mxSn2BDaVdDXpGp4xEfHpvOw1wLcjYoqkBcB3gH8Cfgu8GbhZ0tMRsXdrPwGzoauswrNGzlNBGtGzcq+d3YEdI+KZPOzMIRHxfI5M3J0LB6R81oci4mhJl5AiFr8svIcLgLkR8TUASe8Ejo+I6bnA1bNWXu/kvN6RwN4R8XT1gpImAZMANt544/58BmaWtftQ6/cR8Ux+LuDrksYBrwEbAuvnefMiolK4ZpCuZq44C7ikUnT6aAlweW8WdEjUrHna/a3Wi4XnhwMjgV1zkXqSZVmqWoHQijtJCfN64c5XWf59FpdbFBFL+tNwM+u/dheeohHA3yLiFUl7A5v0cr2fAv/DslFCq80HRktaSdJGwNgG23oBWLsPbTazfmh7SLTgAuA3kqYDs0jjnfdKRHwnh0LPl3R41ew7gHmkIOhc6t90DNKh1O8kPeGTy2at45BoPzgkap3MIVEzsxpceMysdC48ZlY6Fx4zK13TC0+NQOhJNZYZn6MLSHpvrWXqLd+C9o6SNLcV2zaz2lrxdXrDQGi1iLgauLrHBQegBSFUMxuA0g61JL07j/p5O/C+wvSJks7Mz8+T9D1JdyqNo14c9eF1kq6U9ICkH+dgKDnoWdnWoZLOK2zrO/nWGP8paaSk30u6V9JZkh7NmTCAYZLOkXS/pOslrdHij8NsSGtF4Vmj6lBrQo4znAO8B9gLeFOD9TcgpcsPJA2+VzEW+DdgB2BzCsWrga2AfSPi34CvAjdFxC7AlUAx6bkl8IOI2I40KOD7qzfkkUTNmqcVhWdhRIwu/FwMbE0Kej4c6YrFXzZY/6qIeC0iHmBZSBTSyKGP5GzVhaTi1JNLC1msPYGLACLiWuAfheUahVDJ65wdEWMiYszIkSN7sWszq6fMb7V6e4l0MRCqButHjenVQdFiCFXU1yiEamZNVlbheZB0Y67N8+sP9WMbYyVtms/tTCDdGAzgSUnb5OmH1F+d20ljpyNpP9INx8ysDco4x/PNiFhEuonWb/PJ5Uf7sd27SOd85pJCn1fm6ScB15CGQn6iwfqnAPtJuhfYPy/7Qj/aYWYDNGRCopJWA5ZExKt5bPUf9eVr/yKHRK2TdUNIdCidy9gYuCQfkr0MHN3m9pgNWUOm8ETEw8DO7W6HmTmrZWZt4MJjZqVz4TGz0vWq8Eh6k6SLJP0pZ6X+R9JWrW6cpMmSjm/1fsysXD0WHkkiXTMzJSI2j4htgS+xfJzBzKzXetPj2Rt4JSJ+XJmQc023Szpd0lxJcyRNgKX3zrm1TpJ8P0l35YT4pZKG5+nzJZ2Sp8+RtHVh/9tKmpLT6sdVJkq6StKMnCifVJh+lKSH8jrnFJLvIyVdLmla/tkjT58s6dxa+zCz1uhN4dmeFJys9j5gNLATsC9wuqQN8rwVkuT5FhRfJqXFdwGmA58vbO/pPP1HQPHwamvS2OZjga9KWiVPPzIidgXGAMdJWlfSm4GvAG8D3pXXrfhv4IyI2I2UPv9JL/axlNPpZs0zkOt49gQuzOnvJyXdAuwGPE9OkgNIqiTJFwHbAnekozdWJcUgKq7IjzNY/pYXv42IxcBiSX8jHeI9Rio2lWzWRqRbW7wJuKUyLLKkS0m3xoBUHLfN+4Z0f5+1e9jHUh7C2Kx5elN47gcOrTG9Udq7VpJcpLHS6wVEKwnx6nT4CslxSeNJhWT3iHhJ0hRSMr1Rm1bKyy9c7k2kQuR0ulmJenOodROwmqSlEQNJu5HuZzNB0jBJI4FxwNS8SK0k+d3AHpK2yNtYcwDfjI0A/pGLztakQyvy/t8h6fVKwxkXb+h1PfDpwnvoV07LzAaux8KTb9x1CPCu/HX6/cBk4FfAbOA+UnE6MSL+mldbIUkeEU8BE4ELJc0mFaLiOZi+uJbU85kNnJa3RUQ8DnwduAe4AXgAeC6vcxwwRtJsSQ8Ax/Rz32Y2QE1Pp+fDoOMj4sCmbrj3+x8eEQtyj+dK4NyIuLKn9frC6XTrZN2QTh+MVy5PljSLZb2tq9rcHjOr0vSTqBExBZjS7O32Yf++0tmsww3GHo+ZdTgXHjMrXVsLT1/Cp/JQw2aDRtsKj8OnZkNXO3s8fQqfFklaXdLP8vyZkvbO09eUdEm+VudiSfdIGpPnLZD0NUn3Sbpb0vp5es3wqJm1TjsLT3/CpxWfAoiIHUhjdP1caZjkT5KuaN6RdGHhroV11gLujoidgFtZdrP3RuHRpRwSNWueTjy5vDR8GhFPApXwafUy5wNExIOkcbq2YvlhiueSrqyueJk0/hYsP0zxvsCZ+dqfq1k+PLqUhzA2a552hiH7Ez7taZlG674Syy7TLgZBa4ZHzax12tnj6U/4tOJW4PC8zlakMbP+yPLDFG9Luh9QTxweNStZ2wpPP8OnFT8EhkmaA1wMTMz30/khMDKHR7+Qt/McjTk8alayQTWEsaRhwCoRsUjS5sCNwFYR8XIz9+OQqHWybgiJDrYbXq0J3JxvXSrg2GYXHTMbuEFVeCLiBdI9mM2sg3Xi1+lmNsi58JhZ6Vx4zKx0XVN4+pJkz8sv6MU272xuK82sN7qi8LQqyR4Rb29G+8ysb7qi8FA/yT5T0o1aNvTxQdUrShpeb5lKr0jSBkrDLs/Kqfi9ynhTZkNVt3ydXi/Jvgg4JCKeVxoi+W5JV8fyV0X2ZpkPA9dFxNfyRYhrVu9IaXz2SQAbb7xxk96W2dDULYWnHgFflzQOeA3YkHT49dc+LjMNODdfeHhV7k0tx0MYmzVPtxxq3c/y99apOBwYCewaEaOBJ0lDGfdpmYi4lRRGfRw4X9IRzW2+mRV1S+Gpl2TfBPhbRLyS70K4SY11R/S0jKTKds4Bfgrs0oo3YWZJVxxqRURIOgT4rqSTSOdt5pPS7N+TNB2YBTxYY/ULgN/0sMx44ARJrwALAPd4zFpoUKXTy+J0unWybkind8uhlpkNIi48ZlY6Fx4zK50Lj5mVri2FR9KSHE+4L0cZ2pKZ6k2Q1Myar11fpy/MF/Mh6Z+AbwDvKC4gaVhELGlH48ystTrhUOt1pCFtkDRe0s2SfgXMkTRK0tzKgpKOlzQ5Pz8u3x5jtqSL8rR35J7UrDy08dp5+gl5eOLZkk6pboBDomblalePZ408cufqwAbAOwvzxgLbR8Q8SaMabOMkYNOIWCxpnTzteOBTEXGHpOHAIkn7AVvm7Qq4WtK4HJOocEjUrETt6vEsjIjREbE18G7gF/meOwBTI2JeL7YxG7hA0keAV/O0O4DvSDoOWCciXgX2yz8zgXuBrUmFqGga8PHcm9oh3zR+OR7C2Kx52n6oFRF3AeuRgpwALxZmv8rybSyGOw8AfkAKj86QtHJEfBP4F2AN0u0vtib1cr6RC93oiNgiIn5a1QaHRM1K1PbCk4vDMODvNWY/CbxR0rqSVgMOzOusBGwUETcDJwLrAMMlbR4RcyLiP4HppN7NdcCR+dALSRtKemNVGxwSNStRu8/xQOqRfCwiliw72kpyovxU4B5gHssCnsOAX0oakdc/IyKelXRaTqAvAR4AfpfPAW0D3JW3vwD4CPC3wq7G45CoWWkcEu0Hh0StkzkkamZWgwuPmZXOhcfMSufCY2alc+Exs9J1feHplKS7mfVeV9zsvQc9Jt3ryTENRcRrLWyfmVXp+h5PlaVJd6idSs+J9z9I+iEpu7WRpAWSvpZ7TXdLGtCY7GbW2GAoPGvkQ60HgZ8ApwFUpdJHA7vm0UQB3gL8IiJ2johHgbWAuyNiJ+BW4OjqnUiaJGm6pOlPPfVU69+V2SA2GApPvaR7o1T6oxFxd2EbLwPX5OczgFHVO3E63ax5BsM5nqUi4i5JlaR7JZV+VnGZfI+fF6tWfSWWZUeWMMg+F7NOMxh6PEtVJd17TKWbWXsMhr/sNZPuwPV1Uum+j7NZmzmd3g9Op1snczrdzKwGFx4zK50Lj5mVzoXHzErXNYWnGWFQSfPzdT5IurP5rTSz3uimr9MHHAYtTosIp9jN2qRrejxV+hUGLW5A0oL8OF7SFEmXSXpQ0gWFwQXNrAW6qcdTc9jjekMUA/9HCoN+PCI+mZett+2dge2Av5BGI90DuL24gIcwNmueburxNCMMWs/UiHgs35dnFg6JmrVUN/V4lhpAGLSexYXnDomatVg39XiWchjUrLt10192h0HNBgmHRPvBIVHrZA6JmpnV4MJjZqVz4TGz0rnwmFnpWlJ4yh7dM8cerul5STPrBK36Ot2je5pZXWUcai0NdEoaLunG3AuaI+mgPL3W6J7vzsvdJ+nGvNxaks7NgdCZlfWLJL1B0lU5MHq3pB3z9Ml53SmSHpF0XGGdqyTNkHR/zmSZWQu1qsdTM9AJLAIOiYjnoaiK+AAACvNJREFUc+ThbklX53lLA52SRgLnAOMiYp6kN+Rl/h24KSKOlLQOMFXSDVX7PgWYGREHS3on8AvSSKKQclx7A2sDf5T0o4h4BTgyIp6RtAYwTdLlEfH34kYdEjVrnjIOtXYnBTq3J11x/PWcHn8N2BCojFNeDHS+Dbg1IuYBRMQzefp+wHslHZ9frw5UV4E9gffn9W6StK6kEXnebyNiMbBY0t/yvh8DjpN0SF5mI1LIdLnCExFnA2dDuoCwPx+Kmf3/9u4/xo6qDOP49+GngggChUApFCuGFEwKRaXW0m6CiKgBpUgLxFaJDQYRUP+oCUEQJUUjEEQgGGvBALUiVVBpabDtVizQbSy0RZbW0shatFQEiWix9PWPc1Zmb+/+uHd3597dPp9kc+fOnpl7zib7ZmbuPHOSQY9MVAQ6z8qv4yPiv5I2k4oHdA10Cqj2zy3g3Iho77JSOryizS7dyK+7hEElTQFOByZExOuSlhX6ZGaDYNCv8VQEOg8Etuai0wIc081mK4HJko7N++g81VoMXNb5oC5JJ1XZthW4MP9+CrAtIv7ZQxcPBP6Ri87xpKMtMxtEg32NBwqBTkn3AA9JaiM99+bZahtHxEv5msoDkvYAtgIfAa4DbgaezsVnM/CJis2vAX4s6WngdWBGL31dBFyS27cDfXl+j5n1g0OidXBI1JqZQ6JmZlW48JhZ6YbSg8CaxurVUMs8FD6bNevKRzxmVjoXHjMrXVMWHkmfkhT5vprOLNcFhd/PlHRr43poZv3RlIUHmE6aUG9afj8auKDb1mY2pDRd4cnT1EwELuatwjMHmJSf8XNlXnekpEWSNkj6TmH7MyStzMn2nxWmvdks6dpCMr7zaGpy3u+anHg/oLzRmu2emq7wAOcAiyLiOeBlSScDs4EVeSbRm3K7ccD5wPuA8yWNypmwq4DTI+JkoA34SmHf2/L624HOoOnXgEtzqHUS8O9qnZI0S1Jbuuv6pQEdsNnuphkLz3Rgfl6en99X82hEvBoR/wGeIeW+TgXGAo/lyMYMuubBHsivq3lrmuLHgBvz83kOiogd1T6sOIVxyrmaWb2a6j4eSYeQnt1zoqQghUsD+E2V5tWmHRawJCK6K1bbK9oTEXMk/ZqUnH9c0ukRUTVDZmYDo9mOeKYCd0fEMRExOiJGAc+Tnt3Tl2svjwMTJb0HQNJ+kt7b0waSxkTE2oi4gXRqdnz/hmBmvWm2wjMdWFix7ueki8w78mNQr9x1syQiXgJmAvfltPnj9F5IrpC0TtJTpOs7D9fbeTPrG6fT6+B0ujUzp9PNzKpw4TGz0jXVt1pDRa3pdLOBNByujviIx8xK19DCUxkG7aXtFZL2K6NfZja4Gn3EUxkG7ckVwIAWHkk+1TRrgIYVnmphUElT8hTD90t6VtI9Sr4MHAkslbRU0mck3Zi3uVzSprw8RtLv8vJ4Scvz1MSLJR2R1y+TdL2k5cDlkuZJukPSCknPSaqctcLMBlgjj3iqhUEBTiId3YwF3g1MjIhbgC1AS0S0kObOmpTbTwL+LmkkaRbRFZL2Br4PTI2I8cBc4NuFzz4oIiZHxPfy+9HAZODjwB2SdpnQzyFRs4HTyMLTXRj0yYjoiIidpLm3RlduGBF/Bd6RH2ExCrgXOI1UhFaQ5mE/EViSw6JXAUcVdvHTil0uiIidEbEB2ESVu50dEjUbOA25xtFLGLRa+LOalcDnSJPwrQA+D0wAvkqaT319REzoZtt/Vbyv/IJyGHxhada8GnXE010Y9MM9bPMaXYOiraRn6bQCfwBagO0R8SqpGI2QNAFA0t6STuhh3+dJ2kPSGNLpXXsPbc2snxpVeLoLg/b0eNM7gYclLc3vV5BOs1oj4k3gBdI3ZETEG6TidkMOf64BPtTDvtuB5aSA6CX5GT9mNkh2+5CopHnAryLi/r5vc0qkJ2iYla+3f9mhEBL1fSx1GD8eHE43q99uX3giYmaj+2C2u2n0nctmthty4TGz0rnwmFnp6i48tSTL+7CvcZLOGoD9zJM0tb/7MbPB1Z8jnlqS5b0ZR5peps+cLDcbuuoqPD0ky5dLWpBT3nMkXSjpyTxl8Jjc7rzOWR0ktUraB/gmaTbQNZLOl7S/pLmSVuVphc/O287M0xI/BDySk+u3Snomz411WKGPV+ft10m6U0rPDMzp9JvyZ/9R0vslPaA0FfK3+vG3NLO+ioiaf4CLgB/l5d8DJwNTgFeAI4B9gb8A1+Y2lwM35+W1wMi8fFB+nQncWtj/9cBFnW2A54D9c7sO4OD8u08DS0hZryPz50/Nvzu4sL+fAJ/My8uAGwr92lLocwdwSDdjnkW6a7Dt6KOPDrNmBbRFHf/XZf7Ue6rVXbJ8VUS8GBHbgT8Bj+T1a+k6ZfA8SV/IBaOaM4DZOVm+DHgbKfgJaabQl/PyacB9EfFmRGwBflvYR4ukJyStJQVSi1mtBwv9Wl/o8yZSDGMXUUinjxjhdLpZf9R8naSGZPnOwvudvDVl8CWSPkh69s0aSeOqfQxwbkR0CWvm7XpLlpOfp3MbcEpEvCDpGlLx6lTsV2Wffe3IbJDVc8RTT7L8/5SmDH4iIq4GtpGOMCqT54uBywrXZU7qZnetwDRJe+YnDLbk9Z1FZlu+HuVvusyaSD2Fp55kedF388XmdaTC8RSwFBjbeXEZuA7YG3g6t7uum30tBDaQTpluJyXMiYhXgB/m9b8AVvWxb2ZWgt0+nV4PT2FszWwopNN957KZlc6Fx8xK58JjZqVz4TGz0rnwmFnpXHjMrHQuPGZWOhceMyudC4+Zlc6Fx8xK58JjZqVz4TGz0rnwmFnpXHjMrHQuPGZWOhceMyudC4+Zlc6Fx8xK58JjZqVz4TGz0rnwmFnpXHjMrHQuPGZWOhceMyudJ/Srg6TXgPZeGw5th5KmmB7OhusYj4mIEY3uRE/2anQHhqj2Zp+psb8ktXmMNlh8qmVmpXPhMbPSufDU585Gd6AEHqMNGl9cNrPS+YjHzErnwmNmpXPhqZGkMyW1S9ooaXaj+1MvSZslrZW0RlJbXnewpCWSNuTXdxXafz2PuV3SRxvX8+5Jmitpq6R1hXU1j0nS+Py32SjpFkkqeyzDnQtPDSTtCfwA+BgwFpguaWxje9UvLRExrnAvy2zg0Yg4Dng0vyePcRpwAnAmcFv+WzSbeaT+FdUzptuBWcBx+adyn9ZPLjy1+QCwMSI2RcQbwHzg7Ab3aSCdDdyVl+8Czimsnx8R2yPieWAj6W/RVCKiFXi5YnVNY5J0BPDOiFgZ6ZuXuwvb2ABx4anNSOCFwvuOvG4oCuARSaslzcrrDo+IFwHy62F5/VAed61jGpmXK9fbAHJkojbVzvWH6v0IEyNii6TDgCWSnu2h7XAad6fuxjQcx9p0fMRTmw5gVOH9UcCWBvWlXyJiS37dCiwknTr9LZ9qkF+35uZDedy1jqkjL1eutwHkwlObVcBxko6VtA/p4uSDDe5TzSTtL+mAzmXgDGAdaSwzcrMZwC/z8oPANEn7SjqWdMH1yXJ7XbeaxpRPx16TdGr+NuuzhW1sgPhUqwYRsUPSl4DFwJ7A3IhY3+Bu1eNwYGH+lngv4N6IWCRpFbBA0sXAn4HzACJivaQFwDPADuDSiHizMV3vnqT7gCnAoZI6gG8Ac6h9TF8kfUP2duDh/GMDyJEJMyudT7XMrHQuPGZWOhceMyudC4+Zlc6Fx8xK58JjZqVz4TGz0v0PIr+aHTEHqYIAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Show some visualization via tables and/or plots to present your solution.\n", "#TODO (optional): add more visualization, print your solution etc\n", "\n", "# Make a vector with all the shipments\n", "totalshipments = [sum(value(m.shipments[i,j]) for j in cities) for i in cities]\n", "\n", "plt.figure(figsize=(3,10))\n", "\n", "# make a horizontal bar plots\n", "plt.barh(range(cities.size),totalshipments, color=\"blue\",align=\"center\") \n", "\n", "# Use city names on the vertical axis\n", "plt.yticks(range(cities.size),cities,rotation=0)\n", "\n", "# Set title and show plot\n", "plt.title('Total shipments per DC',fontsize=20)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "## Additional visualization (not part of assignment)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Additional packages required. You might need to install these first to make the rest of the code work." ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "# mapping package\n", "import folium\n", "# geocoder\n", "from geopy.geocoders import Nominatim\n", "# limit geocode requests\n", "from geopy.extra.rate_limiter import RateLimiter\n", "# progress bar\n", "from tqdm import tqdm\n", "# namedtuple objects\n", "from collections import namedtuple\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First we find coordinates for each location. This process is called *geocoding*." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "100%|██████████████████████████████████████████████████████████████████████████████████| 37/37 [00:38<00:00, 1.05s/it]\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
locationlatlong
name
Amsterdam(Amsterdam, Noord-Holland, Nederland, (52.3727...52.3727604.893604
Antwerp(Antwerpen, Vlaanderen, België - Belgique - Be...51.2211104.399708
Athens, Greece(Αθήνα, Δήμος Αθηναίων, Περιφερειακή Ενότητα Κ...37.98394123.728305
Barcelona(Barcelona, Barcelonès, Barcelona, Catalunya, ...41.3828942.177432
Berlin(Berlin, 10117, Deutschland, (52.5170365, 13.3...52.51703713.388860
\n", "
" ], "text/plain": [ " location lat \\\n", "name \n", "Amsterdam (Amsterdam, Noord-Holland, Nederland, (52.3727... 52.372760 \n", "Antwerp (Antwerpen, Vlaanderen, België - Belgique - Be... 51.221110 \n", "Athens, Greece (Αθήνα, Δήμος Αθηναίων, Περιφερειακή Ενότητα Κ... 37.983941 \n", "Barcelona (Barcelona, Barcelonès, Barcelona, Catalunya, ... 41.382894 \n", "Berlin (Berlin, 10117, Deutschland, (52.5170365, 13.3... 52.517037 \n", "\n", " long \n", "name \n", "Amsterdam 4.893604 \n", "Antwerp 4.399708 \n", "Athens, Greece 23.728305 \n", "Barcelona 2.177432 \n", "Berlin 13.388860 " ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# an open source geocoder\n", "geolocator = Nominatim(user_agent=\"specify_your_app_name_here\")\n", "# delay as to not overload the server. Note that a commercial one (like GoogleMaps API can handle much more)\n", "geocode = RateLimiter(geolocator.geocode, min_delay_seconds=1)\n", "\n", "# display progress bar\n", "tqdm.pandas()\n", "\n", "# Make a new dataframe\n", "df_cities = pd.DataFrame({\"name\": cities})\n", "df_cities.loc[\n", " df_cities.name == \"Athens\", \"name\"\n", "] = \"Athens, Greece\" # Our open source geocoder is not perfect.. Apparently there is an athens in the US it will use\n", "\n", "# start geocoding\n", "df_cities[\"location\"] = df_cities[\"name\"].progress_apply(geocode)\n", "df_cities[\"lat\"] = df_cities[\"location\"].apply(lambda x: x.latitude)\n", "df_cities[\"long\"] = df_cities[\"location\"].apply(lambda x: x.longitude)\n", "\n", "# set index and display new dataframe with coordinates\n", "df_cities.set_index(\"name\", inplace=True)\n", "df_cities.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A lof of the code below is taken from this blog: https://medium.com/@bobhaffner/folium-lines-with-arrows-25a0fe88e4e" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "def get_arrows(locations, color=\"blue\", size=3, n_arrows=3):\n", "\n", " \"\"\"\n", " Get a list of correctly placed and rotated \n", " arrows/markers to be plotted\n", " \n", " Parameters\n", " locations : list of lists of lat lons that represent the \n", " start and end of the line. \n", " eg [[41.1132, -96.1993],[41.3810, -95.8021]]\n", " arrow_color : default is 'blue'\n", " size : default is 6\n", " n_arrows : number of arrows to create. default is 3\n", " Return\n", " list of arrows/markers\n", " \"\"\"\n", "\n", " Point = namedtuple(\"Point\", field_names=[\"lat\", \"lon\"])\n", "\n", " # creating point from our Point named tuple\n", " p1 = Point(locations[0][0], locations[0][1])\n", " p2 = Point(locations[1][0], locations[1][1])\n", "\n", " # getting the rotation needed for our marker.\n", " # Subtracting 90 to account for the marker's orientation\n", " # of due East(get_bearing returns North)\n", " rotation = get_bearing(p1, p2) - 90\n", "\n", " # get an evenly space list of lats and lons for our arrows\n", " # note that I'm discarding the first and last for aesthetics\n", " # as I'm using markers to denote the start and end\n", " arrow_lats = np.linspace(p1.lat, p2.lat, n_arrows + 2)[1 : n_arrows + 1]\n", " arrow_lons = np.linspace(p1.lon, p2.lon, n_arrows + 2)[1 : n_arrows + 1]\n", "\n", " arrows = []\n", "\n", " # creating each \"arrow\" and appending them to our arrows list\n", " for points in zip(arrow_lats, arrow_lons):\n", " arrows.append(\n", " folium.RegularPolygonMarker(\n", " location=points,\n", " fill_color=color,\n", " number_of_sides=3,\n", " radius=size,\n", " rotation=rotation,\n", " ).add_to(some_map)\n", " )\n", " return arrows\n", "\n", "def get_bearing(p1, p2):\n", "\n", " \"\"\"\n", " Returns compass bearing from p1 to p2\n", " \n", " Parameters\n", " p1 : namedtuple with lat lon\n", " p2 : namedtuple with lat lon\n", " \n", " Return\n", " compass bearing of type float\n", " \n", " Notes\n", " Based on https://gist.github.com/jeromer/2005586\n", " \"\"\"\n", "\n", " long_diff = np.radians(p2.lon - p1.lon)\n", "\n", " lat1 = np.radians(p1.lat)\n", " lat2 = np.radians(p2.lat)\n", "\n", " x = np.sin(long_diff) * np.cos(lat2)\n", " y = np.cos(lat1) * np.sin(lat2) - (np.sin(lat1) * np.cos(lat2) * np.cos(long_diff))\n", " bearing = np.degrees(np.arctan2(x, y))\n", "\n", " # adjusting for compass bearing\n", " if bearing < 0:\n", " return bearing + 360\n", " return bearing" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Produce the map" ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "execution_count": 49, "metadata": {}, "output_type": "execute_result" } ], "source": [ "some_map = folium.Map(tiles='OpenStreetMap')\n", "cent_coord = [df_cities.loc[\"Amsterdam\"][\"lat\"], df_cities.loc[\"Amsterdam\"][\"long\"]]\n", "bounds = []\n", "\n", "for i in cities:\n", " for j in cities:\n", " if value(m.shipments[i, j]) != 0:\n", " origin = i\n", " destination = j\n", " if origin == 'Athens':\n", " origin = 'Athens, Greece'\n", " if destination == 'Athens':\n", " destination = 'Athens, Greece'\n", " origin_coord = [df_cities.loc[origin][\"lat\"], df_cities.loc[origin][\"long\"]]\n", " destination_coord = [df_cities.loc[destination][\"lat\"], df_cities.loc[destination][\"long\"]]\n", " bounds.append([origin_coord])\n", " folium.Circle(\n", " radius=25000,\n", " location=destination_coord,\n", " popup=i,\n", " color=\"#3186cc\",\n", " fill=True,\n", " fill_color=\"#3186cc\",\n", " ).add_to(some_map)\n", " folium.Circle(\n", " radius=80000,\n", " location=origin_coord,\n", " popup=j,\n", " color=\"crimson\",\n", " fill=False,\n", " fill_color=None,\n", " ).add_to(some_map)\n", " # folium.Marker(location=da_coord, popup=da, icon=folium.Icon(color='red')).add_to(some_map)\n", " folium.PolyLine(locations=[origin_coord,destination_coord], color=\"grey\").add_to(\n", " some_map\n", " )\n", " arrows = get_arrows(locations=[origin_coord, destination_coord], n_arrows=1)\n", " for arrow in arrows:\n", " arrow.add_to(some_map)\n", "some_map.fit_bounds(\n", " bounds,\n", " padding_top_left=None,\n", " padding_bottom_right=None,\n", " padding=None,\n", " max_zoom=None,\n", ")\n", "\n", "folium.TileLayer('cartodbdark_matter').add_to(some_map)\n", "folium.TileLayer('cartodbpositron').add_to(some_map)\n", "# other mapping code (e.g. lines, markers etc.)\n", "folium.LayerControl().add_to(some_map)\n", "\n", "some_map" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.1" } }, "nbformat": 4, "nbformat_minor": 2 }