"
],
"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",
"
location
\n",
"
lat
\n",
"
long
\n",
"
\n",
"
\n",
"
name
\n",
"
\n",
"
\n",
"
\n",
"
\n",
" \n",
" \n",
"
\n",
"
Amsterdam
\n",
"
(Amsterdam, Noord-Holland, Nederland, (52.3727...
\n",
"
52.372760
\n",
"
4.893604
\n",
"
\n",
"
\n",
"
Antwerp
\n",
"
(Antwerpen, Vlaanderen, België - Belgique - Be...
\n",
"
51.221110
\n",
"
4.399708
\n",
"
\n",
"
\n",
"
Athens, Greece
\n",
"
(Αθήνα, Δήμος Αθηναίων, Περιφερειακή Ενότητα Κ...
\n",
"
37.983941
\n",
"
23.728305
\n",
"
\n",
"
\n",
"
Barcelona
\n",
"
(Barcelona, Barcelonès, Barcelona, Catalunya, ...
\n",
"
41.382894
\n",
"
2.177432
\n",
"
\n",
"
\n",
"
Berlin
\n",
"
(Berlin, 10117, Deutschland, (52.5170365, 13.3...
\n",
"
52.517037
\n",
"
13.388860
\n",
"
\n",
" \n",
"
\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": [
"