{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"toc": true
},
"source": [
"Table of Contents
\n",
""
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"## Package to install once by running this cell before starting\n",
"import Pkg;\n",
"Pkg.add(\"MetaGraphs\") \n",
"Pkg.add(\"LightGraphs\")\n",
"Pkg.add(\"GraphPlot\")\n",
"Pkg.add(\"JuMP\")\n",
"Pkg.add(\"Ipopt\")\n",
"Pkg.add(\"PyPlot\")"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# Run this cell and start to work while its compiling (will take some time)\n",
"using LightGraphs, MetaGraphs # Graph Libraries\n",
"using JuMP, Ipopt # Optim Libraries\n",
"using GraphPlot, PyPlot # Plot library"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# TP Braess Paradox and Price of Anarchy\n",
"\n",
"## 1 Recall on the user-equilibrium and the social-optimum\n",
"\n",
"Consider a (finite) directed graph $G=(V,E)$. We consider \n",
"$K$ origin-destination vertex pair $\\{o^k,d^k\\}_{k\\in [1,K]}$. \n",
"\n",
"Let denotes by\n",
"+ $r^k$ the intensity of the flow of user entering in $o^k$ and exiting in \n",
"$d^k$;\n",
"+ $\\mathcal{P}_k$ the set of all simple (i.e. \n",
"without cycle) path form $o^k$ to $d^k$, and by $\\mathcal{P}=\\bigcup_{k=1}^K \n",
"\\mathcal{P}_k$;\n",
"+ $f_p$ the flux of user taking path $p \\in \\mathcal{P}$;\n",
"+ $f = \\{f_p\\}_{p\\in\\mathcal{P}}$ the vector of path-flux;\n",
"+ $x_e$ the flux of user taking the edge $e \\in E$;\n",
"+ $x=\\{x_e\\}_{e\\in E}$ the vector of edge-flux;\n",
"+ $\\ell_e : \\mathbb{R} \\to \\mathbb{R}^+$ the cost incurred by a given user \n",
"to take edge\n",
"$e$;\n",
"+ $L_e(x_e) := \\int_0^{x_e}\\ell_e(u)du$.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Minimizing the total cost of the system is an optimization problem that reads\n",
"\n",
"\n",
"\\begin{align}\n",
"(P^{SO}) \\quad \\min_{x,f} \\quad & \\sum_{e\\in E} x_e \\ell_e(x_e) \\\\\n",
" s.t. \\quad & r_k = \\sum_{p\\in\\mathcal{P}_k}f_p & k = 1..K \n",
"\\label{cst:flux-spread}\\\\\n",
" & x_e = \\sum_{p \\ni e} f_p & e \\in E \\label{cst:xa-def}\\\\\n",
" & f_p \\geq 0 & p \\in \\mathcal{P}\n",
"\\end{align}\n",
"\n",
"Where the first constraint ensure that the flux going from $o^k$ \n",
"to $d^k$ is spread among the different possible paths \n",
"and the second constraint is the definition of $x_e$ as the sum of the \n",
"users taking the different path containing edge $e$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Recall that if $\\ell_e$ is non-decreasing we can guarantee that $f$ is a user \n",
"equilibrium (cf course) if and only if it is a solution of\n",
"\n",
"\\begin{align}\n",
"(P^{UE}) \\quad \\min_{x,f} \\quad & \\sum_{e\\in E}L_e(x_e) \\\\\n",
" s.t. \\quad & r_k = \\sum_{p\\in\\mathcal{P}_k}f_p & k \\in 1..K\\\\\n",
" & x_e = \\sum_{p \\ni e} f_p & e \\in E \\\\\n",
" & f_p \\geq 0 & p \\in \\mathcal{P}\n",
"\\end{align}\n",
"where the only difference with $(P^{SO})$ is the \n",
"objective function."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 1 \n",
"Prove that if for all edge $e$ the cost $\\ell_e$ is constant then social optimum and \n",
"user-equilibrium are equivalent."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 2 Reformulating the problem\n",
"\n",
"A complete graph $(V,E)$ is a graph such that for all couple of vertices\n",
"$(v_1,v_2)$ there exists an edge with origin $v_1$ and destination $v_2$.\n",
"\n",
"### Question 2\n",
"\n",
"Consider a complete directed graph of $n$ nodes with $1$ origin and $1$ \n",
"destination. What is the dimension of the vector $f$ ? of $x$ ? So how many \n",
"variables and constraints (except positivity constraints) is there for the user \n",
"equilibrium (or the social optimum) problem ?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 3\n",
"\n",
"Rewrite both $(P^{(UE)})$ and $(P^{(SO)})$ using only vector $f$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 4\n",
"\n",
"Assuming that there is only one origin and destination (i.e. $K=1$)\n",
" Rewrite both problems using only vector $x$. We can use, for any node $i$ the notation $in(i)$ for the set of edges with destination $i$, and $out(i)$ for the set of edges with origin $i$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 5\n",
"\n",
"Let $N$ be a matrix with $|V|$ line and $|E|$ column, such that $N_{od}$ \n",
"equal to $-1$ if the origin of $e$ is $o$, $+1$ if the destination of $e$ is \n",
"$d$ and $0$ elsewhere. Thus each column of $N$ corresponds to an edge and \n",
"indicate its origin with a $-1$ and its destination with a $+1$.\n",
"\n",
"If $K=1$, what is the interpretation of each coordinate of the vector $Nx$ ? \n",
"Deduce a more condensed presentation of both problems using only vector $x$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setting up the optimization problem in Julia"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We now set up the graph we want to study, which is the classical Braess example."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{4, 4} directed Int64 metagraph with Float64 weights defined by :weight (default weight 1.0)"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# A utility function to add an edge o->d with metadata (a,b,idx) to the graph\n",
"import MetaGraphs.add_edge!\n",
"\n",
"function add_edge!(G::AbstractMetaGraph,o,d,a,b,idx)\n",
" add_edge!(G, o, d) # adding edge o -> d\n",
" set_props!(G, o, d,Dict(:a=>a, :b => b, # the cost of edge o -> d is a x + b = x \n",
" :idx => idx)) # idx stands for the index corresponding to the edge-flow in the x vector\n",
"end\n",
"\n",
"\n",
"# Constructing a Directed Graph (with Meta-data) with 4 nodes\n",
"\n",
"\n",
"function create_graph()\n",
" G = MetaDiGraph(4)\n",
"\n",
" # setting the input / output flux\n",
" set_props!(G,1,Dict(:io=>-1,:x=>0,:y=>0)) # incoming flux of intensity 1 \n",
" set_props!(G,2,Dict(:io=>0,:x=>1,:y=>1))\n",
" set_props!(G,3,Dict(:io=>0,:x=>1,:y=>-1))\n",
" set_props!(G,4,Dict(:io=>1,:x=>2,:y=>0)) # outgoing flux of intensity 1\n",
"\n",
" # adding edges with attached metadata\n",
" add_edge!(G, 1, 2, 1., 0., 1) # edge 1->2 with cost function 1.*x+0. and idx 1\n",
" add_edge!(G, 1, 3, 0., 1., 2)\n",
" add_edge!(G, 2, 4, 0., 1., 3)\n",
" add_edge!(G, 3, 4, 1., 0., 4)\n",
"\n",
" return G\n",
"end\n",
"\n",
"G = create_graph()"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGgCAYAAABMn6ZGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xl4VOX5xvHvZE9IyEJMwhICYd+FsEMUESGAlLhCRRQKKi1UEWwpYlGjgruoFBRlUVSKUlFahILIKgEDEtklyBKWJEAC2deZ+f2RH6fGgLJkcibJ/bmuuep7cnLmHtsyD+9zzvta7Ha7HREREZFqysXsACIiIiKOpGJHREREqjUVOyIiIlKtqdgRERGRak3FjoiIiFRrKnZERESkWlOxIyIiItWaih0RERGp1lTsiIiISLWmYkdERESqNRU7IiIiUq25mR3ADDabjdOnT+Pn54fFYjE7joiIiFwBu91OdnY29erVw8XlyudramSxc/r0acLDw82OISIiItfgxIkTNGjQ4IrPr5HFjp+fH1D6L6t27dompxEREZErkZWVRXh4uPE9fqVqZLFzsXVVu3ZtFTsiIiJVzNXegqIblEVERKRaU7EjIiIi1ZqKHREREanWVOyIiIhItaZiR0RERKo1FTsiIiJSranYEZFqaebMmXTp0gU/Pz9CQkKIjY3lxx9/NDuWiJhAxY6IVEsbN25k/PjxbNu2jbVr11JSUkL//v3Jzc01O5qIVDKL3W63mx2ismVlZeHv709mZqYWFRSpIc6ePUtISAgbN27kpptuMjuOiFyDa/3+1syOiFQrBcVWzmYXUlBsLXM8MzMTgKCgIDNiiYiJauR2ESJS/SQcy+D9zUdYuz8Nmx1cLHBb61Aeio4kKiKQSZMm0bt3b9q2bWt2VBGpZA6d2dm0aRNDhgyhXr16WCwWvvjii9/8nY0bNxIVFYWXlxeRkZG888475c6ZM2cOjRs3xsvLi6ioKDZv3uyI+CJSRSzedpx734nn6wNnsP1/Y95mh68PnOGed+K57e4H2b17N0uWLDE3qIiYwqHFTm5uLh06dGD27NlXdP7Ro0cZNGgQ0dHR7Nq1iyeffJJHH32Uf/3rX8Y5S5cuZeLEiUybNo1du3YRHR3NwIEDSU5OdtTHEBEnlnAsg+lf7MUOWG1lb0G02uykr32HDWtX8ebi5TRo0MCckCJiqkq7QdlisbB8+XJiY2Mve86UKVNYsWIFBw4cMI6NGzeOH374gfj4eAC6detGp06dmDt3rnFOq1atiI2NZebMmVeURTcoi1QfjyzewdcHzpQrdOx2O+e/foe8Q/HUG/Eig3p34p37o0xKKSIV4Vq/v53qnp34+Hj69+9f5tiAAQOYP38+xcXF2O12du7cyd/+9rcy5/Tv35+tW7de9rqFhYUUFhYa46ysrIoNLiKmKCi2GvfoAJTkZJC17TOKzyXj6h9G3o9bCLnzKexu3ny1/QDHbg4jNDgIb29vc4OLSKVyqmInNTWV0NDQMsdCQ0MpKSnh3Llz2O12rFbrJc9JTU297HVnzpzJs88+65DMImKe7IISrMVF5CVtI2fvOgqO7gK77f9/+gMAaUumGuc3ng0LFy5k1KhRlR9WREzjVMUOlLa7fu5il81isZT551+e88tjPzd16lQmTZpkjLOysggPD6+oyCJSyex2O9u2beP9+Qs48dES7IW5eNRrgcXdExd3L3xa30x2QukDEX5RQwi89WFcXSzsj4vBy93V5PQiUtmcqtgJCwsrN0Nz5swZ3NzcqFOnDna7HVdX10ue88vZnp/z9PTE09PTIZlFpPIkJyezePFi5i9cxNGfDuMZEIJfp8H4tulL4akDpK96k6AhT+DTtBu+bfuSsvBRsnf+m+yd/2b0vI0qdERqKKdaVLBHjx6sXbu2zLE1a9bQuXNn3N3d8fDwICoqqtw5a9eupWfPnpUZVUQqSW5uLh9++CF9bulLo0aNeDruec54hRMy7HlCH36fwJsewOLuSca696jV9lZ8mnYDwCMkkoaTP8fFyxeAhQ/fzMaNG838KCJiEofO7OTk5HD48GFjfPToURITEwkKCqJhw4ZMnTqVU6dO8eGHHwKlT17Nnj2bSZMm8dBDDxEfH8/8+fPLrI0xadIkRo4cSefOnenRowfz5s0jOTmZcePGOfKjiEglstlsbNq0iUWLPuDTzz4jPy8Xn4h2BMU8ik+LXrh4+hjn2u12Mla/jYuHN8G3PVzmOm4enjR87J90ydzIZ++8Qp8+ffjDH/7A+++//6utbxGpXhz66PmGDRu45ZZbyh1/8MEHWbRoEaNGjeLYsWNs2LDB+NnGjRt5/PHH2bdvH/Xq1WPKlCnlCpk5c+bw8ssvk5KSQtu2bXnjjTeuaq8bPXou4ryKi4uJbNqMk8nH8Qyqi1frvtRqcwvuAWGXPD87cTUZ/53Nm4s+Zb9rJGv2pRorKPdvE8bY3o3p3CiIAwcO0Lp1a+P3zp8/T0BAQGV9LBGpANf6/a2NQFXsiDiV8+fP06JlKzIuZFLnjqfwbnTjZc8tyTxD2sLxPHj/fbz//vtA6ePo2QUl+Hm5lbtHp6ioiGbNmhmLkP73v/8tt9yFiDgvbQQqItVCYGAge3b/wM29e3Hm079zYfNH2G3WcufZ7TbOr36TkOA6vP7668ZxL3dXbvDzvOTNyB4eHhw/fpxZs2YBpet4DRs2jBr4dz6RGkXFjog4ndDQUNauXcPT06eTufWfJL8ylJLsc2XOydm1irxjP/DBooVXPUP72GOPkZSUBMCnn36Ki4sL6enpFZZfRJyLih0RcUqHDh0qsxjomQ8eI/+nBACKL6SSuXEh48aNo1+/ftd0/aZNm1JUVESbNm0ACA4OZsWKFdcfXEScjoodEXE6L730Eq1atQJg1KhRnD17lltv6smZZc9yfv0CLqyaRd2wUF555ZXreh93d3f27t3LvHnzABg6dCiDBg3CZrP9xm+KSFWiG5R1g7KI0ygoKKB+/fpkZGQAsH79evr06QOUPo7+xhtvMOVvf8NaUlLmZxUhOTmZiIgIY5ySkkJY2KWfABMRc+gGZRGp0nbv3o23t7dR6GRmZpYpZlxcXJg8eTIJ333H5s2bK7TQAWjYsCHFxcV061a6KGHdunVZunRphb6HiJhDxY6ImG769Ol06NABgAkTJmCz2S77t7aOHTvSu3dvh+Rwc3Nj27ZtfPTRRwAMHz6cm266Cau1/NNgIlJ1qI2lNpaIafLy8ggICKC4uBiArVu30qNHD5NTlTp9+jT169c3xidOnKBBgwYmJhIRtbFEpEpJSEigVq1aRqGTk5PjNIUOQL169bBarcbTXuHh4SxatMjcUCJyTVTsiEilstvtTJo0ia5duwIwZcoU7HY7tWrVMjlZeS4uLqxdu5Zly5YBMHr0aDp27EhJSYnJyUTkaqiNpTaWSKXJycnBz8/PGO/cuZNOnTqZmOjKnTlzhtDQUGN85MgRGjdubGIikZpHbSwRcWpbtmwxCp1atWqRl5dXZQodgJCQEKxWK3fccQcAkZGRzJkzx+RUInIlVOyIiEPZ7XYeeeQRoqOjAYiLiyMnJwdvb2+Tk109FxcXPv/8c1auXAnA+PHjadasGUVFRSYnE5FfozaW2lgiDpOZmUlAQIAx3rt3r7E9Q1WXkZFBnTp1jPHBgwdp0aKFiYlEqj+1sUTEqaxbt84odOrWrUtBQUG1KXQAgoKCsNlsjBw5EoCWLVvy6quvmpxKRC5FxY6IVCi73c6IESOMR7Zfe+01Tp8+jaenp8nJKp7FYuHDDz9k3bp1APzlL38hLCyMwsJCk5OJyM+5mR1ARKqP9PR0goODjfGPP/5I8+bNTUxUOfr27UtmZib+/v6kpaXh5eXF7t27adeundnRRATN7IhIBVm5cqVR6DRv3pyioqIaUehcVLt2bWw2G+PGjQOgffv2PPvssyanEhFQsSMi18lmszF06FBuv/12AObOncuPP/6Iu7u7yckqn8ViYe7cuWzZsgWAZ555Bh8fH/Lz801OJlKzqY0lItfslwvtHT16lEaNGpkXyEn06tWL7Oxs/Pz8yM/Px8fHhx07dhAVFWV2NJEaSTM7InJNli1bZhQ6UVFRFBcXq9D5GV9fX+x2O0888QQAnTt35q9//avJqURqJq2zo3V2RK6K1Wqlf//+fPPNNwAsWrSIBx980ORUzm3Hjh106dLFGOfk5DjlXmAizu5av7/VxhKRK3b69Gnq169vjE+ePFlmLJfWuXNn8vLyCAoKoqCgAF9fX7799lt69uxpdjSRGkFtLBG5IosXLzYKmz59+lBSUqJC5yp4e3uTn5/P008/DZTe1/OnP/2JGji5LlLp1MZSG0vkV5WUlNCrVy++++47AJYuXcq9995rcqqqbc+ePbRv394Y688ikSuj7SJEpMIdP34cd3d3o9BJTU1VoVMB2rVrR35+PjfccAMA/v7+rF+/3uRUItWXih0RuaR58+YZT1fdfvvtWK3WMo+Zy/Xx8vLizJkzvPTSS0DpKswPPvig2loiDqA2lqaORcooLi6mffv2HDx4EIAVK1YwZMgQk1NVbz/++CMtW7Y0xhkZGQQGBpqYSMQ5qY0lItctKSkJDw8Po9A5d+6cCp1K0KJFCwoLC2nSpAlQuqP6qlWrTE4lUn2o2BERAGbNmmXsZTV8+HBsNht16tQxOVXN4eHhweHDh5k9ezYAgwYN4s4771RbS6QCqI2lNpbUcBdnFE6dOgXAmjVruO2220xOVbMdOXLEmOWB0m05Lt7MLFKTOXUba86cOTRu3BgvLy+ioqLYvHnzZc/t06cPFoul3Gvw4MHGOaNGjSr38+7du1fGRxGpVvbv34+Xl5dR6Jw/f16FjhOIjIykuLiYDh06ABASEsLy5ctNTiVSdTm82Fm6dCkTJ05k2rRp7Nq1i+joaAYOHEhycvIlz//8889JSUkxXnv37sXV1ZV77rmnzHkxMTFlzvvqq68c/VFEqpUXXniBNm3aADB27FhsNhsBAQEmp5KL3NzcSExMZMGCBQDceeedDBgwAJvNZnIykarH4W2sbt260alTJ+bOnWsca9WqFbGxscycOfM3f3/WrFlMnz6dlJQUYy+ZUaNGceHCBb744otryqQ2ltRk+fn5hIWFkZWVBcCmTZuIjo42OZX8mhMnTtCwYUNjfPr0aerWrWtiIhFzOGUbq6ioiJ07d9K/f/8yx/v378/WrVuv6Brz589n+PDh5TbN27BhAyEhITRv3pyHHnqIM2fOXPYahYWFZGVllXmJ1ESJiYn4+PgY/x/IyspSoVMFhIeHU1JSQu/evQGoV68en3zyicmpRKoOhxY7586du+RCZKGhoaSmpv7m73/33Xfs3buXsWPHljk+cOBAPv74Y7755htee+01EhIS6Nu3L4WFhZe8zsyZM/H39zde4eHh1/6hRKqoJ598ko4dOwIwceJEbDYbfn5+JqeSK+Xq6srmzZtZsmQJACNGjKBnz55YrVaTk4k4P4e2sS7ukLx161Z69OhhHH/hhRdYvHixsZbH5TzyyCNs3bqVPXv2/Op5KSkpRERE8M9//pM777yz3M8LCwvLFEJZWVmEh4erjSU1Qm5uLn5+fsYjzNu3b6dr164mp5LrkZqaWqaNdfz48TJtLpHqyinbWMHBwbi6upabxTlz5sxvLjufl5fHP//5z3KzOpdSt25dIiIiSEpKuuTPPT09qV27dpmXSE2wfft2fH19sdvtuLq6kpubq0KnGggLC8NqtRITEwNAREQE8+fPNzmViPNyaLHj4eFBVFQUa9euLXN87dq19OzZ81d/99NPP6WwsJD777//N98nPT2dEydO6IY9kf9nt9t57LHHjCUZnnzySUpKSvDx8TE5mVQUFxcXVq1aZTySPnbsWNq1a0dJSYnJyUScj5uj32DSpEmMHDmSzp0706NHD+bNm0dycjLjxo0D4IEHHqB+/frlnsyaP38+sbGx5VZwzcnJ4ZlnnuGuu+6ibt26HDt2jCeffJLg4GDuuOMOR38cEaeXnZ1dZvZy165d3HjjjSYmEkeKjY3l7Nmz3HDDDezduxd3d3cOHz5cZlFCkZrO4evsDBs2jFmzZhEXF8eNN97Ipk2b+Oqrr4iIiAAgOTmZlJSUMr9z6NAhtmzZwpgxY8pdz9XVlT179jB06FCaN2/Ogw8+SPPmzYmPj9fNllLjbdy40Sh0AgICyM/PV6FTAwQHB2Oz2bj77rsBaNq0KW+//bbJqUSch7aL0P07Ug3Y7XbGjh1rLED3wgsv8OSTT5qcSsywevVqBg4cCJTey3Po0CE8PDxMTiVSMa71+9vhbSwRcawLFy4QGBhojPfv30+rVq1MTCRmiomJISMjg6CgII4fP46npycHDhygZcuWZkcTMY12PRepwtasWWMUOuHh4RQWFqrQEQIDA7HZbIwePRooXbX+xRdfNDmViHlU7IhUQXa7neHDhzNgwACgdFuV5ORktSvEYLFYWLBgAevXrwdg6tSp1KlTh4KCApOTiVQ+tbFEqpj09HSCg4ONcVJSEk2bNjUxkTizPn36kJWVRe3atcnIyMDb25vExERjR3WRmkAzOyJVyIoVK4xCp02bNhQVFanQkd/k5+eHzWZjwoQJANx4441Mnz7d5FQilUfFjkgVYLPZGDx4MEOHDgVg3rx5xpoqIlfCYrHw9ttvG5swP/fcc3h4eJCXl2dyMhHHUxtLxMlpHySpSD169CAnJwdfX1+Ki4upVasW3333HV26dDE7mojDaGZHxIktXbrUKHS6detGcXGxCh25brVq1cJutzNlyhQAunbtyqRJk6iBy65JDaFFBbWooDghq9VK37592bRpEwAfffQRI0aMMDmVVEfff/89UVFRxjg7OxtfX18TE4lcnhYVFKkmTp48SXh4uDE+deoU9erVMzGRVGedOnUiLy+PkJAQcnJy8PPzY/PmzfTu3dvsaCIVRm0sESeyaNEio9Dp168fVqtVhY44nLe3N9nZ2cTFxQEQHR3NI488oraWVBtqY6mNJU6gpKSErl27smvXLgCWLVvGXXfdZXIqqYn27dtH27ZtjfGFCxfw9/c3MZHI/1zr97dmdkRMdvToUdzd3Y1CJy0tTYWOmKZNmzYUFBQYN8YHBASwbt06k1OJXB8VOyImmjNnDpGRkQDExsZitVoJCQkxOZXUdJ6enpw+fZrXXnsNKG2pjhgxQm0tqbLUxlIbS0xQVFREmzZtOHz4MAArV65k0KBBJqcSKS8pKYnmzZsb4/T0dIKCgkxMJDWZ2lgiVcShQ4fw9PQ0Cp309HQVOuK0mjVrRlFRES1atACgTp06/Oc//zE5lcjVUbEjUolee+0140tj5MiR2Gw2/S1ZnJ67uzsHDx5k7ty5AAwZMoTf/e532Gw2k5OJXBkVO1LjbNq0iSFDhlCvXj0sFgtffPGFw9+zsLCQsLAwnnjiCQDWrVvHhx9+iMVicfh7i1SUcePGcfToUQD+/e9/4+rqypkzZ0xOBXPnzqV9+/bUrl2b2rVr06NHD1atWmV2LHEiKnakxsnNzaVDhw7Mnj27Ut5v7969eHl5kZaWBpQ+ytu3b99KeW+RitaoUSOKi4vp3LkzAKGhoSxbtszUTA0aNODFF19kx44d7Nixg759+zJ06FD27dtnai5xHip2pMYZOHAgzz//PHfeeafD3ysuLo527doBpX8rttlsWrNEqjw3NzcSEhL44IMPALjnnnvo27cvVqvVlDxDhgxh0KBBNG/enObNm/PCCy/g6+vLtm3bTMkjzkfbRYg4QH5+PsHBweTl5QGwZcsWevXqZXIqkYr1wAMPcOutt9KgQQPWr1+Pm5sbJ0+epH79+qZlslqtfPbZZ+Tm5tKjRw/Tcohz0cyOSAX7/vvv8fHxMQqd7OxsFTpSbdWvX5+SkhL69OkDlLaUFi9eXOk59uzZg6+vL56enowbN47ly5fTunXrSs8hzknFjkgFmjJlirGD9OTJk7Hb7dpBWqo9V1dX1q9fz9KlS4HSGZ+uXbtSUlJSaRlatGhBYmIi27Zt449//CMPPvgg+/fvr7T3F+emRQW1qGCNZrFYWL58ObGxsdd1ndzc3DJFTUJCgnEDp0hNkpaWRlhYmDE+duwYERERlZ6jX79+NGnShHfffbfS31scR4sKipgkPj7eKHQ8PT3Jy8tToSM1VmhoKFarlSFDhgClT2+ZUXDY7XYKCwsr/X3FOanYkRonJyeHxMREEhMTgdKNOBMTE0lOTr6q69jtdsaPH0/Pnj0BmD59OgUFBXh7e1d4ZpGqxMXFhRUrVrBixQqg9EnEVq1aUVxc7JD3e/LJJ9m8eTPHjh1jz549TJs2jQ0bNjBixAiHvJ9UPWpjqY1V42zYsIFbbrml3PEHH3yQRYsWXdE1Lv5v6KLdu3cbj5iLyP+kp6cTHBxsjJOSkmjatGmFvseYMWNYt24dKSkp+Pv70759e6ZMmcJtt91Woe8j5rvW728VOyp25CqtX7/eWBQwODiYEydO4OXlZXIqEedlt9u57777+Oc//wnAG2+8wcSJE01OJVWR7tkRcTC73c6oUaOMQuell17i7NmzKnREfoPFYmHJkiWsWbMGgMcff5wGDRronhqpNFpUUOQKnD9/vsyGnQcPHjQ29BSRK3Pbbbdx4cIFAgICOHXqFF5eXuzbt0/r4YjDaWZH5DesXr3aKHQiIyMpLCxUoSNyjfz9/bHZbIwdOxaANm3a8MILL5icSqq7Sil25syZQ+PGjfHy8iIqKorNmzdf9txFixZhsVjKvQoKCq75miLXwm63c9dddzFw4EAAZs+ezU8//YSHh4fJyUSqNovFwnvvvcemTZsAeOqpp/D39y/357xIRXF4sbN06VImTpzItGnT2LVrF9HR0QwcOPBXH/OtXbs2KSkpZV4/vy/iWq4pcjXOnj2Li4sLn3/+OQA//fQT48ePNzmVSPUSHR1NVlYWUHrjqbe3N7t27TI5lVRHDi92Xn/9dcaMGcPYsWNp1aoVs2bNIjw8nLlz5172dywWC2FhYWVe13PNwsJCsrKyyrxELmf58uWEhIQA0KFDB4qLi4mMjDQ5lUj15Ofnh81mM57O6tSpE08++aTJqaS6cWixU1RUxM6dO+nfv3+Z4/3792fr1q2X/b2cnBwiIiJo0KABt99+e5lK/1quOXPmTPz9/Y1XeHj4dXwqqa5sNhsDBgzgzjvvBGDBggUkJibi5qb7+EUcyWKx8MYbb7B9+3ag9M9sFxcXcnNzTU4m1YVDi51z585htVoJDQ0tczw0NJTU1NRL/k7Lli1ZtGgRK1asYMmSJXh5edGrVy+SkpKu+ZpTp04lMzPTeJ04caICPp1UJykpKbi6uhqPxiYnJzN69GiTU4nULF27diU3NxdXV1djE92LBZDI9aiUG5QtFkuZsd1uL3fsou7du3P//ffToUMHoqOj+fTTT2nevDlvv/32NV/T09OT2rVrl3mJXPTJJ59Qr149AHr37k1JSYlm/0RM4uPjQ0lJCdOmTQNKvxMee+wxauD6t1KBHFrsBAcH4+rqWm7G5cyZM+VmZi7HxcWFLl26GDM7FXFNEQCr1UqvXr2M/XM++eQTNm/ejKurq8nJROT555839q976623cHFxITs72+RUUlU5tNjx8PAgKiqKtWvXljm+du1aY/PE32K320lMTKRu3boVdk2REydO4ObmZtznlZKSwu9//3uTU4nIz3Xo0IH8/HwCAgKA0id1Lz6uLnI1HN7GmjRpEu+//z4LFizgwIEDPP744yQnJzNu3DgAHnjgAaZOnWqc/+yzz/Lf//6XI0eOkJiYyJgxY0hMTDTOv5Jrivya+fPn07BhQwBiYmKwWq3lnvgTEefg5eXF+fPnmTFjBgA333wzY8aMUVtLrorDHzMZNmwY6enpxMXFkZKSQtu2bfnqq6+IiIgASm8EdXH5X8114cIFHn74YVJTU/H396djx45s2rSJrl27XvE1RS6lpKSETp06sWfPHqD0EfPY2FiTU4nIlZg6dSqxsbG0bt2aBQsWsGDBAs6fP2/M+oj8Gu16rpuVa4QjR47QpEkTY3z27FmCg4NNTCQi16KoqIimTZsaT9X+97//LbcUiVRf2vVc5DLefvtto9C5++67sdlsKnREqigPDw+Sk5OZNWsWAAMGDGDYsGFqa8mv0syOZnaqraKiIlq0aMGxY8cAWLVqFTExMeaGEpEKc/jwYZo1a2aMz507R506dUxMJI6mmR2Rnzl48CCenp5GoZORkaFCR6Saadq0KUVFRbRp0wYoXZpkxYoVJqcSZ6RiR6qdl156iVatWgEwatQobDYbgYGBJqcSEUdwd3dn7969zJs3D4ChQ4cyaNAgbDabycnEmaiNpTZWtVFQUECDBg1IT08HYP369fTp08fcUCJSaZKTk8s8lZuamqrFZqsZtbGkRtu9ezfe3t5GoZOVlaVCR6SGadiwIcXFxXTr1g2AsLAwli5danIqcQYqdqTKe/rpp+nQoQMAEyZMwGaz4efnZ3IqETGDm5sb27Zt46OPPgJg+PDh3HTTTVitVpOTiZnUxlIbq8rKy8sjICCA4uJiALZu3UqPHj1MTiUizuL06dPUr1/fGJ84cYIGDRqYmEiul9pYUqMkJCRQq1Yto9DJyclRoSMiZdSrVw+r1Uq/fv0ACA8PZ9GiReaGElOo2JEqZ/Lkycb2IVOmTMFut1OrVi2TU4mIM3JxcWHt2rUsW7YMgNGjR9OpUydKSkpMTiaVSW0stbGqjJycnDL34uzcuZNOnTqZmEhEqpIzZ86UeTrryJEjNG7c2MREcrXUxpJqbcuWLUah4+vrS15engodEbkqISEhWK1WYwPgyMhI5syZY3IqqQwqdsSp2e12xo0bR3R0NABxcXFkZ2fj7e1tcjIRqYpcXFxYvnw5K1euBGD8+PE0a9bMuP9Pqic3swOIXE5mZiYBAQHGeO/evcay8CIi12PQoEGkp6dTp04dDh8+jIeHBwcPHqRFixZmRxMH0MyOOKV169YZhU7dunUpKChQoSMiFSooKAibzcbIkSMBaNmyJa+++qrJqcQRVOyvW9qPAAAgAElEQVSIU7Hb7YwYMcJ4VPTVV1/l9OnTeHp6mpxMRKoji8XChx9+yLp16wD4y1/+QlhYGIWFhSYnk4qkNpY4jYyMDOrUqWOMDx06RLNmzUxMJCI1Rd++fcnMzMTf35+0tDS8vLzYs2cPbdu2NTuaVADN7IhTWLlypVHotGjRgqKiIhU6IlKpateujc1mY9y4cQC0a9eOuLg4k1NJRVCxI6ay2WzExsZy++23AzB37lwOHjyIu7u7yclEpCayWCzMnTuXLVu2AKV77/n4+JCfn29yMrkeamOJaX65wNfRo0dp1KiReYFERP5fr169yM7Oxs/Pj/z8fHx8fNixYwdRUVFmR5NroJkdMcW//vUvo9CJioqiuLhYhY6IOBVfX1/sdjuTJ08GoHPnzkyZMsXkVHIttF2EtouoVFarlf79+/PNN98A8MEHH/DAAw+YnEpE5Nft2LGDLl26GOOcnBztyWeCa/3+VhtLKs3p06epX7++MT558mSZsYiIs+rcuTN5eXkEBgZSWFiIr68v3377LT179jQ7mlwBtbGkUixevNgobPr06UNJSYkKHRGpUry9vSkoKGD69OlA6X0948ePpwY2SKoctbHUxnKokpISevbsSUJCAgBLly7l3nvvNTmViMj12bNnD+3btzfG+j6pHNr1XJzO8ePHcXd3Nwqd1NRUFToiUi20a9eO/Px8goODAfD392f9+vUmp5LLUbEjDjFv3jzj6aohQ4ZgtVrLPGYuIlLVeXl5cfbsWV566SWgdBXmBx98UG0tJ6Q2lqYdK1RxcTEdOnTgwIEDAKxYsYIhQ4aYnEpExLF+/PFHWrZsaYwzMjIIDAw0MVH1pDaWmO7w4cN4eHgYhc65c+dU6IhIjdCiRQsKCwuJjIwESndUX716tcmp5CIVO1IhZs2aZexlNXz4cGw2W5lNPUVEqjsPDw9++uknZs+eDcDAgQO566671NZyApVS7MyZM4fGjRvj5eVFVFQUmzdvvuy57733HtHR0QQGBhIYGEi/fv347rvvypwzatQoLBZLmVf37t0d/THkEoqKiggPD+fxxx8HYM2aNSxZsgSLxWJyMhERc4wfP56ffvoJgM8//xwXFxfOnj1rcqqazeHFztKlS5k4cSLTpk1j165dREdHM3DgQJKTky95/oYNG/j973/P+vXriY+Pp2HDhvTv359Tp06VOS8mJoaUlBTj9dVXXzn6o8gv7N+/H09PT06ePAnAhQsXuO2220xOJSJivsjISOMeRoCQkBCWL19ucqqay+HFzuuvv86YMWMYO3YsrVq1YtasWYSHhzN37txLnv/xxx/zpz/9iRtvvJGWLVvy3nvvYbPZWLduXZnzPD09CQsLM15BQUGO/ijyMy+88AJt2rQBYOzYsdhsNvz9/U1OJSLiPNzc3EhMTGTBggUA3HnnnQwYMACbzWZysprHocVOUVERO3fupH///mWO9+/fn61bt17RNfLy8iguLi5XzGzYsIGQkBCaN2/OQw89xJkzZy57jcLCQrKyssq85NoUFBQQEBDAU089BcCmTZt477331LYSEbmM0aNHG92MNWvW4OrqSkpKismpahaHFjvnzp275PoqoaGhpKamXtE1/va3v1G/fn369etnHBs4cCAff/wx33zzDa+99hoJCQn07duXwsLCS15j5syZ+Pv7G6/w8PBr/1A1WGJiIt7e3mRmZgKljwBGR0ebnEpExPmFh4dTUlJC7969AahXrx6ffPKJyalqjkq5QfmXf+u32+1XNBPw8ssvs2TJEj7//HO8vLyM48OGDWPw4MG0bduWIUOGsGrVKg4dOsTKlSsveZ2pU6eSmZlpvE6cOHF9H6gGmjZtGh07dgRg4sSJ2Gw2/Pz8TE4lIlJ1uLq6snnzZqPIGTFiBD179sRqtZqcrPpz6K7nwcHBuLq6lpvFOXPmzG+upvvqq68yY8YMvv766zL7j1xK3bp1iYiIICkp6ZI/9/T0xNPT8+rCC1DaRvT19TUendy+fTtdu3Y1OZWISNX1+9//nltuuYW6desSHx+Pm5sbycnJ6jo4kENndjw8PIiKimLt2rVljq9du5aePXte9vdeeeUVnnvuOVavXk3nzp1/833S09M5ceIEdevWve7M8j/bt2+nVq1a2O12XF1dyc3NVaEjIlIBwsLCsFqtxMTEANCwYUPmz59vcqrqy+FtrEmTJvH++++zYMECDhw4wOOPP05ycjLjxo0D4IEHHmDq1KnG+S+//DJPPfUUCxYsoFGjRqSmppKamkpOTg4AOTk5PPHEE8THx3Ps2DE2bNjAkCFDCA4O5o477nD0x6kR7HY7jz32mLF20bRp0ygpKcHHx8fkZCIi1YeLiwurVq0yHkkfO3Ys7dq1o6SkxORk1Y9D21hQen9Neno6cXFxpKSk0LZtW7766isiIiIASE5OxsXlfzXXnDlzKCoq4u677y5znaeffppnnnkGV1dX9uzZw4cffsiFCxeoW7cut9xyC0uXLtU9JBUgOzu7zH4jiYmJxjoRIiJS8WJjYzl79iw33HADe/fuxd3dncOHD9OkSROzo1Ub2ghUG4EaNm3axM033wxAQEAAKSkpZW4MFxERx7Hb7dx7770sW7YMgLfeeos///nPJqdyLtoIVK6Z3W5n7NixRqEzY8YMzp8/r0JHRKQSWSwWPvvsM1atWgXAo48+SuPGjSkqKjI5WdXn8DaWOLcLFy4QGBhojPfv30+rVq1MTCQiUrPFxMSQkZFBUFAQx44dw9PTkwMHDtCyZUuzo1VZmtmpwdasWWMUOg0bNqSwsFCFjoiIEwgMDMRmszFq1CgAWrVqxYsvvmhuqCpMxU4NZLfbGTZsGAMGDABg1qxZHD9+HA8PD5OTiYjIRRaLhYULF7J+/XqgdIHcOnXqUFBQYHKyqkdtrBomPT2d4OBgY5yUlETTpk1NTCQiIr+mT58+ZGVlUbt2bTIyMvD29uaHH374zQV35X80s1ODrFixwih02rRpQ1FRkQodEZEqwM/PD5vNxoQJEwDo0KED06dPNzlV1aFipwaw2WwMHjyYoUOHAjBv3jxjLQcREakaLBYLb7/9Nlu3bgXgueeew8PDg7y8PJOTOT+1saq5tLQ0wsLCjPHx48dp2LChiYlEROR69OjRg5ycHHx9fSkuLqZWrVp89913dOnSxexoTkszO9XY0qVLjUKne/fulJSUqNAREakGLu5bOGXKFAC6du3K5MmTTU7lvLSCcjVcQdlqtdK3b182bdoEwEcffcSIESNMTiUiIo7w/fffExUVZYyzs7Px9fU1MZHjXOv3t9pY1czJkycJDw83xqdOnaJevXomJhIREUfq1KkTeXl5hISEkJOTg5+fH5s3b6Z3795mR3MaamNVIx988IFR6PTr1w+r1apCR0SkBvD29iY7O5u4uDgAoqOjeeSRR6iBzZtLUhurGrSxSkpK6Nq1K7t27QJg2bJl3HXXXSanEhERM+zbt4+2bdsa4wsXLuDv729iooqjjUBrqKNHj+Lu7m4UOmlpaSp0RERqsDZt2lBQUEDdunUBCAgIYN26dSanMpeKnSps7ty5REZGAhAbG4vVaiUkJMTkVCIiYjZPT09Onz7Nq6++CpTe2jBixIga29ZSG6sKtrGKi4tp06YNSUlJAKxcuZJBgwaZnEpERJxRUlISzZs3N8bp6ekEBQWZmOjaqY1VQxw6dAgPDw+j0ElPT1ehIyIil9WsWTOKioqMgqdOnTqsXLnS5FSVS8VOFfLaa6/RokULAEaOHInNZquy1bmIiFQed3d3fvzxR+bOnQvA7bffztChQ7HZbCYnqxxqY1WBNlZhYSERERGkpaUBsG7dOvr27WtyKhERqYqOHTtG48aNjXFaWlqVud9Tbaxqau/evXh5eRmFTmZmpgodERG5Zo0aNaK4uNhYdTk0NJRly5aZnMqxVOw4sbi4ONq1awfAH//4R2w2W5WYiRIREefm5ubGjh07+OCDDwC455576Nu3L1ar1eRkjqE2lhMWD/n5+QQHB5OXlwfAt99+S8+ePU1OJSIi1dGpU6do0KBBmbGzrr6vNlY18f333+Pj42MUOtnZ2Sp0RETEYerXr09JSQl9+vQxxosXLy533unTp7l9yO/4w5ixVW69HhU7TuRvf/ub0UN94oknsNvt1XbnWhERcR6urq6sX7+epUuXAvDAAw/QtWtXSkpKAFizZg3t2ndgzfpNLFwwn3fffdfMuFdNbSwnaGPl5uaWKWoSEhLo3LmziYlERKSmSktLIywszBg/9thjvPnmm/hERhE0eBIXNi/GemgT+/buLfNUV2VQG8tJbdq0iSFDhlCvXj0sFgtffPFFmZ/Hx8cbhY6npyd5eXkqdERExDShoaFYrVbju+jNN98k4OZRBN/9NK4+/gT2+QM2Dz8eHDX6itfpmTlzJhaLhYkTJzoy+mWp2HGw3NxcOnTowOzZs8sct9vtjB8/3rgfZ/r06RQUFODt7W1GTBEREcPKlStJOvwTLu6eBMU8in/3u7FYSksGF08f/GMeZfOmjfzjH//4zWslJCQwb9482rdv7+jYl+Vm2jvXEAMHDmTgwIFljl2chrto9+7dxiPmIiIiZikqKmLKlCnMmjWLWs26Ue+Bibh6+5U7zzuiA36dBvOXv05h4MCBNG3a9JLXy8nJYcSIEbz33ns8//zzjo5/WZrZcYCCYitnswspKC6/XsHevXuNQueGG26goKBAhY6IiDiFnTt3MmvWLFx9g6jV8XZcPH0ue27AzaPBJ4CRDzxorM/zy++/8ePHM3jwYPr161cp+S9HMzsVKOFYBu9vPsLa/WnY7OBigdtah/JQdCRREYEA/P3vfwfgpZde4q9//auZcUVERMq48cYbmTFjBvMXLuKnT/+Op/8NeLbqg2/bW3Gv06DMuS4eXgTEPMa2JVN54pkXyWs+oMz3X+OsHzi+PYG9id+b9Gl+ltXsANXF4m3HufedeL4+cAbb/z/fZrPD1wfOcOestbi4/O9f9cGDB1XoiIiI0/H29mbq1Kkk/XiQ+Ph4Rg2/E/v+NZx+fxxnPnqC7F1fYS3IMc73Cm+LX9TvePPFOL7a8r3x/VeUeZaNi14hv+cfWZaYZtKn+Z9KKXbmzJlD48aN8fLyIioqis2bN//q+f/6179o3bo1np6etG7dmuXLl5f5ud1u55lnnqFevXp4e3vTp08f9u3b58iP8KsSjmUw/Yu92AGrreyT/DmHd3DizeHG+NNPPzV2LhcREXFGFouF7t27884773AmLZVPP/2UPh0iubDuXU7/YyTnvnyRvJ8SsNusBNw0Ele/YM7853XsttL2VVHqYax5Fzi9aCIje0bi5ubGxo0beeutt3Bzc6v0bSkcXuwsXbqUiRMnMm3aNHbt2kV0dDQDBw4kOTn5kufHx8czbNgwRo4cyQ8//MDIkSO599572b59u3HOyy+/zOuvv87s2bNJSEggLCyM2267jezsbEd/nEt6f/MRXFwsZY7Z7XbOLp/Bmc+eBqBO/z8C4O7uXun5RERErpWXlxf33HMPX61cyamTJ3n5xZk0dL3A2WXPcnrOKC5s/ojaXe+gKCWJrITSyQmviA7U/cNs6o5+iwZj3ubO5z6mc+fOjBgxgsTERFxdXSv1Mzh8UcFu3brRqVMn5s6daxxr1aoVsbGxzJw5s9z5w4YNIysri1WrVhnHYmJiCAwMZMmSJdjtdurVq8fEiROZMmUKAIWFhYSGhvLSSy/xyCOPlLtmYWEhhYWFxjgrK4vw8PAKWVSwoNhK6+mr+fmETklOBqf+8YAxvuGu6bj5BZOy6FFeeuVV+ve7laCgIBo2bHhd7y0iImIGu93O9oSd9P/Tc+Tu34AtP6v0By6u1B39Nh7BZb/fXCzQYOvLdOrYkVmzZl3z+zrlooJFRUXs3LmT/v37lznev39/tm7desnfiY+PL3f+gAEDjPOPHj1KampqmXM8PT25+eabL3vNmTNn4u/vb7zCw8Ov52OVkV1Qwi86V+Tt31j6Dy6u3HD3M5z9Vxwpix4FYMpfnqBjx45Mnz69wjKIiIhUJovFQpNW7Qjq9zD1Hn4Pvy6xuNQKAJuNguO7y51vs5e/zaMyOfRprHPnzmG1WgkNDS1zPDQ0lNTU1Ev+Tmpq6q+ef/E/L3XO8ePHL3nNqVOnMmnSJGN8cWanIvh5ueFioUzBU6tdP86vnw82K2eXPUPdP8zG44ZGuFhgf1wMXu6VO30nIiJSkex2Oz/siCd91VvkHtyMvSgfsGBxc8e7ccdy57tYYN036037/quUG5QtlvL3s/zy2NWefzXX9PT0pHbt2mVeFcXL3ZXbWofi+rN7dly9/Wj413/je2MMACkLJpC19Z/0bxOmQkdERKqso0eP8uyzz9Iosgm39b2FwuTd1O4SS1DMo4CdgJtH4R5Uv8zvuLpYTP/+c2ixExwcjKura7lZnDNnzpSbmbkoLCzsV8+/uDnZ1VzT0cZGR2L7xfScxWKhzoAJhI54GYDzmz/iw0duJj8/34yIIiIi1yQ7O5uFCxcSfdPNREZG8vzMl0n3a0LofS9S9+F51O5yB5nxS/EMb4tf1O3lft9mszO2d+VuGPpLDi12PDw8iIqKYu3atWWOr1271tgT6pd69OhR7vw1a9YY5zdu3JiwsLAy5xQVFbFx48bLXtPRujQK4rnYtligzAwPQK2GbWj4+GcA5Ofl4uPjw/ffm7/AkoiIyOVYrVa+/vpr7r9/JCGhYfxhzBi+P5lFndsnU/dPHxI8aCJe4W2xWFw4v2EhtrxMQgZPNPbPgtLvQwvwXGxbOjcKMu/DUAkrKE+aNImRI0fSuXNnevTowbx580hOTmbcuHEAPPDAA9SvX994Muuxxx7jpptu4qWXXmLo0KF8+eWXfP3112zZsgXA2DV1xowZNGvWjGbNmjFjxgx8fHy47777HP1xLuv+7hG0DPPj/S1HWbMvtcwKymN7N6bz63aeeOIJXnvtNaKiopgyZQovvviiaXlFRER+KTc3lxkzZrBg0Qeknj6FV3ADvLrcTZ02t+BW+4Zy5+cf3UVO4ir+GvcymY1vvPT3n8mFDlTCo+dQuqjgyy+/TEpKCm3btuWNN97gpptuAqBPnz40atSIRYsWGecvW7aMp556iiNHjtCkSRNeeOEF7rzzTuPndrudZ599lnfffZfz58/TrVs3/vGPf9C2bdsrynOtj65dqYJiK9kFJfh5uZXrUSYkJNC1a1djnJOTQ61atSo8g4iIyNXaunUrvXr1wtU3iOChU/Gs3/Ky98PaCvM4s3AC3Tu24Zt1X+Pi4vKr338V4Vq/vyul2HE2ji52fkteXh6BgYEUFRUBpY/bd+/evdJziIiI/FxJSQlxcXE8//zzeEd0IGjwZFx9Ay95bvrqt7Ad/pb9+/YRERFRKfmccp0duTQfHx8KCwuNtXZ69OjBhAkTqIF1p4iIOBE3Nzfi4uJYu3YttfJSOPmPkWQlfFHuvPwjO8n5YQ2z3nij0gqd66GZHRNmdn5u9+7ddOjQwRg7QyYREanZSkpKaN26NUlJSQD49xiGf+/7sLi4YivIIW3RBG7q0pE1a/77q0vJVDTN7FRR7du3Jz8/n+DgYAD8/f3ZsGGDuaFERKTGOn78OO7u7kah8/e//53s7Z9xbuk0SrLPkbHuPdxthSxYML9SC53roWLHCXh5eXH27Fnj6axbbrmF0aNHq60lIiKVat68eTRq1AiAwYMHY7VaiYuLY8OGDfgVZ5C2YAK5e9fx9ptvVujWS46mNpaTtYwOHjxIq1atjHFGRgaBgZe+OUxERKQiFBcX06FDBw4cOADAl19+ye9+97sy55w7d44Jf36UkNAQ3nzjDVNmdfQ01lVw5mIHShdJbNmyJUePHgVg1apVxMTEmJxKRESqo6SkJJo3b26Mz507R506dUxMdHm6Z6ca8fDw4MiRI7z99tsADBw4kHvuuUdtLRERqVCzZs0yCp3hw4djs9mcttC5HprZccKZnZ+7uLDiRWfPnjVuZhYREbkWhYWFNGnShFOnTgGl2zLddtttJqf6bZrZqaYiIyMpLi6mffv2ANxwww188UX5NQ9ERESuxP79+/Hy8jIKnfPnz1eJQud6qNipAtzc3Pjhhx94//33AbjjjjsYOHAgNpvN5GQiIlKVzJgxgzZt2gAwduxYbDYbAQEBJqdyPLWxnLyN9UsnTpygYcOGxjglJYWwsDATE4mIiLMrKCggLCyMzMxMADZt2kR0dLTJqa6e2lg1RHh4OCUlJfTs2ROAunXrsmTJEpNTiYiIs0pMTMTb29sodLKysqpkoXM9VOxUQa6urnz77bd88sknANx333307t0bq9VqcjIREXEmTz75JB07dgRg4sSJ2Gw2/Pz8TE5V+dTGqmJtrF9KTU2lbt26xjg5OblKrWopIiIVLzc3Fz8/P2PJkm3bttGtWzeTU10/tbFqqLCwMKxWKwMGDACgYcOGLFiwwORUIiJilu3bt+Pr64vdbsfV1ZXc3NxqUehcDxU71YCLiwurV6/m888/B2DMmDF06NCBkpISk5OJiEhlsdvtPPbYY3Tv3h0obWGVlJTg4+NjcjLzqY1VxdtYv3T27FlCQkKM8U8//URkZKSJiURExNGys7PLfJ/t2rWLG2+80cREjqE2lgCliw7abDbuvvtuAJo0acLs2bNNTiUiIo6yadMm44s/ICCA/Pz8alnoXA8VO9WQxWLhs88+Y9WqVQD8+c9/pkmTJhQVFZmcTEREKordbmfMmDHcfPPNALzwwgucP38eLy8vk5M5HzezA4jjxMTEkJGRQVBQEEeOHMHT05MDBw7QsmVLs6OJiMh1uHDhAoGBgcZ4//79tGrVysREzk0zO9VcYGAgNpuNUaNGAdCqVStefvllc0OJiMg1W7NmjVHohIeHU1BQoELnN6jYqQEsFgsLFy5k/fr1AEyZMoUbbriBgoICk5OJiMiVstvtDB8+3FhqZNasWSQnJ+Pp6WlyMuenNlYN0qdPHzIzM/H39+fcuXN4e3vzww8/GDuqi4iIc0pPTyc4ONgYJyUl0bRpUxMTVS2a2alhateujc1mY/z48QB06NCBp59+2uRUIiJyOStWrDAKndatW1NUVKRC5yqp2KmBLBYLs2fPZuvWrQDExcXh5eVFXl6eyclEROQim83G4MGDGTp0KADz5s1j3759uLu7m5ys6lEbqwbr0aMHOTk5+Pr6UlhYSK1atUhISKBz585mRxMRqdHS0tIICwszxseOHSMiIsLERFWbZnZquFq1amG32/nrX/8KQJcuXXjiiSdMTiUiUnMtXbrUKHS6detGcXGxCp3rpO0iqtl2Edfj+++/JyoqyhhnZ2fj6+trYiIRkZrDarXSt29fNm3aBMDixYu5//77TU7lXK71+1ttLDF06tSJvLw8goODycvLw8/Pjy1bttCrVy+zo4mIVGsnT54kPDzcGJ86dYp69eqZmKh6URtLyvD29iY3N5dnn30WgN69ezNu3Dhq4ASgiEilWLRokVHo3HrrrVitVhU6Fcyhxc758+cZOXIk/v7++Pv7M3LkSC5cuHDZ8zMyMvjzn/9MixYt8PHxoWHDhjz66KNkZmaWOc9isZR7vfPOO478KDXO9OnT2bNnDwDvvvsuLi4u5f57EBGRa1dSUkKnTp0YPXo0AMuWLePrr7/GxUXzEBXNoW2s++67j5MnT7J69WoAHn74YUaOHMm///3vS55/+vRpTp8+zauvvkrr1q05fvw448aN4/Tp0yxbtqzMuQsXLiQmJsYY+/v7O+6D1FBt27aloKCARo0akZqaSkBAAOvWraNv375mRxMRqdKOHj1KZGSkMU5LSyMkJMTERNWc3UH2799vB+zbtm0zjsXHx9sB+8GDB6/4Op9++qndw8PDXlxcbBwD7MuXL7/mbJmZmXbAnpmZec3XqGleffVVO2AH7CNGjLDbbDazI4mIVElz5swx/jyNjY21W61WsyNVGdf6/e2wubL4+Hj8/f3p1q2bcax79+74+/sbi9ldiYt3XLu5lZ2EmjBhAsHBwXTp0oV33nkHm8122WsUFhaSlZVV5iVXZ/Lkyfz4448AfPzxx7i4uJCRkWFyKhGRqqO4uJjmzZvzpz/9CYD//Oc/LF++XG2rSuCwf8OpqamXnJILCQkhNTX1iq6Rnp7Oc889xyOPPFLm+HPPPcdnn33G119/zfDhw5k8eTIzZsy47HVmzpxp3Dfk7+9f5o53uXLNmzenqKiIZs2aAVCnTh2++uork1OJiDi/Q4cO4eHhQVJSElD6/TZ48GCTU9UcV13sPPPMM5e8Qfjnrx07dgClNxL/kt1uv+TxX8rKymLw4MG0bt263N5NTz31FD169ODGG29k8uTJxMXF8corr1z2WlOnTiUzM9N4nThx4io/tVzk7u7OoUOHmDNnDgCDBw8mNjb2V2fWRERqstdee40WLVoAcP/992Oz2QgKCjI5Vc1y1TcoT5gwgeHDh//qOY0aNWL37t2kpaWV+9nZs2cJDQ391d/Pzs4mJiYGX19fli9f/pv7gHTv3p2srCzS0tIueW1PT088PT1/9Rpydf74xz8SExNDZGQkX375Ja6urrrBTkTkZwoLC4mIiDC+C7/++mtuvfVWk1PVTFdd7AQHB5fZZv5yevToQWZmJt999x1du3YFYPv27WRmZtKzZ8/L/l5WVhYDBgzA09OTFStW4OXl9ZvvtWvXLry8vAgICLjyDyLXrXHjxhQXF9O9e3d27txJaGgoy5Yt46677jI7moiIqfbu3Uu7du2M8YULF/TUsIkcds9Oq1atiImJ4aGHHmLbtm1s27aNhx56iNtvv92Yzjt16hQtW7bku+++A0pndPr3709ubi7z588nKyuL1NRUUlNTsVqtAPz73//mvffeY+/evfz000+8/2qlQMkAACAASURBVP77TJs2jYcfflizNyZwc3Njx44dLFy4EIC7776bfv36qa0lIjVWXFycUeiMGzcOm82mQsdkDt0bKyMjg0cffZQVK1YA8Lvf/Y7Zs2cbMzDHjh2jcePGrF+/nj59+rBhwwZuueWWS17r6NGjNGrUiNWrVzN16lQOHz6MzWYjMjKSsWPHMn78+HJPbF2O9sZyDC13LiI1WX5+vrHdDqDtdhzgWr+/tRGoip0KZbVaufXWW9m4cSOgjexEpGbQRsqV41q/v/Vwv1QoV1dXNmzYwNKlSwEYOXIk3bt3N9qQIiLVzZQpU4xCZ/LkydjtdhU6TkYzO5rZcZi0tDTCwsKM8bFjx4iIiDAxkYhIxcnNzS1T1CQkJNC5c2cTE1V/mtkRpxMaGorVajUWzmrUqBHvvfeeyalERK5ffHy8Ueh4enqSm5urQseJqdgRh3JxceE///kPX375JVC6GWybNm0oLi42OZmIyNWz2+2MHz/eWEJl+vTpFBQU4OPjY3Iy+TVqY6mNVWnS09PLrNGUlJRE06ZNTUwkInLlLn53XPTDDz/Qvn17ExPVPGpjidOrU6cONpuNYcOGAdCsWTPefPNNk1OJiPy29ev/r707j4riStsA/jRbgwgNitCggGgUxAUVFXAlLqjRGCfjuCQyMqNmEyOjmWjMTCRm4jbZ40SjEnUmbnHBaExwSQDjIEgUoiIS44oKisq+032/PxzrswPIIk11N8/vnD7xVt+qem+K7nq7btW9sVKi4+TkhNLSUiY6RoTJDjUrhUKB7du349ChQwCAiIgIeHh4oKKiQubIiIiqE0IgLCwMw4cPBwCsXLkSOTk59RrdnwxHg6eLIGoKo0aNQm5uLhwdHZGZmQmlUolz586hW7ducodGRAQAyM3N1ZmwMz09HT4+PjJGRI3FKzskGwcHB2i1WsycORMA4Ovri+XLl8scFREREBMTIyU6nTp1Qnl5ORMdI8Zkh2SlUCiwYcMGacTlxYsXw8HBAWVlZTJHRkQtkRACkyZNwtixYwEAn376KS5evAgrKyuZI6PHwW4sMghDhw5FQUEB7O3tkZ+fDxsbG6SkpKB3795yh0ZELUROTg6cnZ2l8sWLF9GpUycZI6Kmwis7ZDDs7Oyg1Woxb948AECfPn3wt7/9TeaoiKgliI6OlhIdPz8/VFZWMtExIRxnh+PsGKSkpCQEBgYCuD/fVkFBAQftIqImp9Vq8dRTT+HgwYMAgA0bNkj3EZLhaez5m91YZJACAgJQXFwMOzs7aDQa2NraIikpCQMGDJA7NCIyEVlZWXBzc5PK165dg7u7u4wRkb6wG4sMVqtWraDRaPDGG28AuJ8AzZs3Dy3wYiQRNbFt27ZJic6gQYNQVVXFRMeEsRuL3VhGITU1FX369JHKBQUFsLOzkzEiIjJGGo0GQ4cORUJCAgBg69atmDZtmsxRUX2xG4tMWu/evVFaWgq1Wi39kR89ehRDhgyROzQiMhKZmZnw8PCQyllZWVCr1TJGRM2F3VhkNKytrZGXl4d3330XwP3H1WfPnm0Q3VrLly9H//79YWdnB2dnZ0ycOBEZGRlyh0Vk0JYvXw6FQoGIiAi97ysqKkpKdEaPHg2NRsNEpwVhskNGZ/HixUhLSwNw/8kJMzMz5OXlyRpTfHw85syZg8TERBw+fBhVVVUICQlBcXGxrHERGark5GSsW7dO75NpVlVVoVevXpg1axaA+4+Yx8TEwMyMp7+WhEebjJKvry/KysrQoUMHAICjoyMOHz4sWzwxMTEICwtD9+7d4efnh40bN+LatWs4efKkbDERGaqioiI8//zzWL9+PRwdHfW2n0uXLsHS0hJnzpwBcH/QwIkTJ+ptf2S4mOyQ0VIqlcjMzMRHH30EAAgJCcG0adMMolsrPz8fAHQmESSi++bMmYNx48Zh5MiRetvH6tWr0blzZwDApEmToNVq4eTkpLf9kWFjskNGb968ebhw4QIAYPv27TAzM8Pdu3dli0cIgfnz52Pw4MHo0aOHbHEQGaLt27fj1KlTepv0t6KiAl5eXpg7dy4A4LvvvsPOnTuhUCj0sj8yDkx2yCQ88cQTqKioQLdu3QAATk5O2L9/vyyxhIeH4/Tp09i2bZss+ycyVJmZmZg3bx6+/PJLWFtbN/n2z58/D6VSiStXrgAA7t27hzFjxjT5fsj4MNkhk2FpaYlz585h3bp1AIAJEyZg/Pjx0Gq1zRbD3LlzsW/fPsTGxkr3ExHRfSdPnsTt27fh7+8PCwsLWFhYID4+Hp988gksLCyg0Wgave2VK1dKP3bCwsKg1Wr1ej8QGRcOKshBBU3S1atX0bFjR6mcnZ0NFxcXve1PCIG5c+ciOjoacXFx6NKli972RWSsCgsLcfXqVZ1lf/rTn+Dj44OFCxc2qtv3wYMKD7quY2NjERwc3BThkgHioIJED/H09ERlZSUGDx6MpKQkqNVq7NixA5MnT9bL/ubMmYOtW7fi66+/hp2dHbKzswEAKpUKNjY2etknkbGxs7OrltDY2tqibdu2jUp0Tp8+DT8/P6nMH7BUG3ZjkcmysLBAYmIi/vOf/wAApkyZguDg4Me6VF6bNWvWID8/H8HBwXB1dZVeO3bsaPJ9ERGwZMkSKdEJDw+HVqtlokO1YjcWPxwtws2bN9G+fXupfP36dZ0yERmHkpISODo6oqKiAgCQkJCAoKAgmaOi5tLY8zev7FCL4ObmBo1GgxEjRgAAOnTogM2bN8scFRE1RHJyMmxtbaVEp6ioiIkO1QuTHWoxzMzMcOTIEezcuRPA/Sc2+vXrh6qqKpkjI6K6LFiwAAMGDAAALFy4EEII2NrayhwVGQt2Y7Ebq0W6ffu2ztNZly9f1nl6i4gMQ1FREezs7KTyyZMn0bdvXxkjIjkZZDdWbm4uQkNDoVKpoFKpEBoaWueEjcHBwVAoFDqvqVOnPvZ2iR7m7OwMjUaDZ555BgDg5eWFNWvWyBwVET3s2LFjUqJja2uLkpISJjrUKHpNdp577jmkpqYiJiYGMTExSE1NRWhoaJ3rzZ49G1lZWdLr888/b5LtEj3MzMwMe/fuxTfffAMAeOWVV+Dt7Y3KykqZIyNq2YQQeOmllzBkyBAAwNKlS1FUVMRhHKjR9DbOTnp6OmJiYpCYmIiAgAAAwPr16xEUFISMjAx4e3vXum6rVq2gVqubfLtENRk3bhzu3r2Ltm3b4pdffoGVlRUyMjLQtWtXuUMjanHy8/Ph4OAglc+cOcM55uix6e3KzvHjx6FSqaSEBAACAwOhUqmQkJDwyHW3bNkCJycndO/eHa+99hoKCwsfa7vl5eUoKCjQeRE9rE2bNtBqtXj++ecBAN7e3nj//fdljoqoZfn++++lRMfV1RVlZWVMdKhJ6C3Zyc7OhrOzc7Xlzs7O0uiyNXn++eexbds2xMXF4e9//zt2796NZ5999rG2u3z5cun+HpVKBXd390a0iEydQqHAl19+iSNHjgAAXnvtNbi5uaG8vFzmyIhMmxAC06dPx8iRIwEA7733Hm7evAmlUilzZGQqGtyNFRkZibfffvuRdZKTkwHcP3n8lhCixuUPzJ49W/p3jx490KVLF/Tr1w+nTp2Sbkxr6HbfeOMNzJ8/XyoXFBQw4aFajRgxAnl5eXBwcEBWVhasra15KZ1IT+7du4e2bdtKZXYhkz40+MpOeHg40tPTH/nq0aMH1Go1bt26VW39nJycBk3I2LdvX1haWuLChQsA0KjtKpVK2Nvb67yIHkWlUkGr1eLFF18EAPTs2RPvvPOOzFERmZYDBw5IiU7Xrl1RUVHBRIf0osHJjpOTE3x8fB75sra2RlBQEPLz83HixAlp3aSkJOTn52PgwIH13l9aWhoqKyvh6uoKAE22XaK6KBQKrF27Fj/++CMA4K233kLr1q1RWloqc2RExk2r1WLixIkYP348AOCzzz5DRkYGLC0tZY6MTJVeBxUcO3Ysbt68KT06/sILL8DT0xP79+8HANy4cQMjRozAv//9bwwYMAAXL17Eli1b8NRTT8HJyQnnzp3DggULYGNjg+TkZJibm9dru3XhoILUUBzYjKhpcEBPehwGOajgli1b0LNnT4SEhCAkJAS9evWSZqAGgMrKSmRkZKCkpAQAYGVlhe+//x6jR4+Gt7c3Xn31VYSEhODIkSNSolOf7RI1tdatW0MIId375e/vj0WLFskcFZFx2b17t5To+Pv7o7KykokONQtOF8ErO9RAycnJ0hw9wP2rPpyjh6h2Go0GISEh+OGHHwAAmzZtwowZM2SOioxRY8/fehtUkMhU9e/fH8XFxXB0dERFRQVat26N48ePIzAwUO7QiAzOzZs30b59e6l8/fp1nTJRc+Cs50SN0KpVK5SXl+Pvf/87gPs3zoeHh6MFXiglqtWXX34pJTbDhg1DVVUVEx2SBbux2I1Fj+n06dPw8/OTyvy7opauqqoKgwYNkp6a3bFjByZPnixzVGQKDPIGZaKWoFevXigtLYWTkxOA+2P0xMXFyRsUkUyuXr0KS0tLKdHJzs5mokOyY7JD1ASsra2Rk5ODFStWAACefPJJhIWFsVuLWpR169ZJT1eNHz8eGo2mQYPIEukLu7HY3UBN7Pz58+jWrZtUvnfvHhwdHWWMiEi/Kisr4efnh/T0dADAvn378PTTT8scFZkidmMRGQgfHx+Ul5dLv3DbtGmDmJgYeYMi0pNff/0VVlZWUqJz584dJjpkcJjsEOmBlZUVLl++jE8++QTA/VG/J02axG4tMikfffQRunTpAgCYMmUKtFqtzqSeRIaC3VjsxiI9u3TpEjp37iyVc3JypJuZiYxRRUUFOnfujOvXrwMADh06hFGjRskcFbUE7MYiMlCdOnVCRUUFevXqBQBo164d9u7dK3NURI1z7tw5KJVKKdHJzc1lokMGj8kOUTOwtLTEzz//jA0bNgAAfve732HMmDHQarUyR0ZUf8uWLUP37t0BALNmzYJWq4WDg4PMURHVjd1Y7MaiZpaZmQkPDw+pfPPmTbi6usoYEdGjlZWVQa1WIz8/HwBw9OhRDBkyROaoqCViNxaRkXB3d5dGmAUANzc3bNu2TeaoiGqWmpoKGxsbKdEpKChgokNGh8kOkQzMzc1x7NgxbNmyBQDw3HPPYfDgwdBoNDJHRvT/3nzzTfTp0wcAEBERAa1WCzs7O5mjImo4dmOxG4tklpWVBTc3N6l87do1uLu7yxgRtXQlJSVo3bq1NFRCUlISBgwYIHNUROzGIjJarq6u0Gg0GD16NADAw8MDX3zxhcxRUUuVlJQEW1tbCCFgbm6O4uJiJjpk9JjsEBkAMzMzxMTEYM+ePQCAmTNnws/PD1VVVTJHRi2FEALz5s1DYGAgAGDx4sWoqqpCq1atZI6M6PGxG4vdWGRgcnJy4OzsLJUvXryITp06yRgRmbrCwkKd78KUlBT07t1bxoiIasZuLCIT0a5dO2i1Wvz+978HAHTu3BmrV6+WOSoyVUePHpVOGg4ODigtLWWiQyaHyQ6RAVIoFNi1axe+/fZbAMDcuXOlkZiJmoIQArNmzcKwYcMAAO+++y5yc3NhbW0tc2RETc9C7gCIqHZjx47FvXv30KZNG1y+fBlKpRLp6enw8fGROzQyYnl5eXB0dJTKaWlp8PX1lTEiIv3ilR0iA+fo6AitVosZM2YAALp164ZVq1bJHBUZq0OHDkmJjru7O8rLy5nokMljskNkBBQKBTZt2oTY2FgAwMKFC+Hs7IyysjKZIyNjIYTA1KlTpSEOPvroI1y7dg1WVlYyR0akf+zGIjIiwcHByM/Ph0qlQk5ODmxsbPDzzz9LM6oT1eTu3btwcnKSyhcuXMATTzwhY0REzYtXdoiMjL29PbRaLebMmQMA8PPzw5IlS2SOigzVvn37pETH19cXFRUVTHSoxWGyQ2SEFAoFVq9ejYSEBADA0qVLoVQqUVJSInNkZCi0Wi3Gjx+PZ555BgCwbt06pKWlwdLSUubIiJofu7GIjFhQUBCKiorQunVrVFRUwNbWFsnJyejXr5/coZGMbt26BbVaLZWvXr0KDw8PGSMikhev7BAZuQfzGP31r38FAPTv3x+vvfaazFGRXL766isp0QkICEBlZSUTHWrxOF0Ep4sgE3Ly5EmdqzqFhYVo3bq1jBFRc9FoNBg+fDiOHj0KAPjPf/6D6dOnyxwVUdNq7Pmb3VhEJsTf3x8lJSVo164diouLYWdnh2PHjmHQoEFyh0Z6dP36dbi7u0vlGzduwM3NTcaIiAwLu7GITIyNjQ2Kiorw9ttvAwAGDx6Ml156CS3wIm6LsHnzZinRGTFiBDQaDRMdot/Qa7KTm5uL0NBQqFQqqFQqhIaGIi8vr9b6V65cgUKhqPG1c+dOqV5N769du1afTSEyOm+99RbOnDkDAPj8889hZmaG/Px8maOiplJVVYW+ffsiLCwMALBr1y4cOXIEZmb8DUv0W3q9Z2fs2LG4fv061q1bBwB44YUX0LFjR+zfv7/G+hqNBjk5OTrL1q1bh1WrViE7O1u690ChUGDjxo0YM2aMVE+lUsHGxqZecfGeHWpJysvL0bFjR2RnZwMAjhw5ghEjRsgcFT2OK1euwMvLSyrfunULzs7OMkZE1Dwae/7W20+A9PR0xMTEYMOGDQgKCkJQUBDWr1+Pb775BhkZGTWuY25uDrVarfOKjo7GlClTqt1k6eDgoFPvUYlOeXk5CgoKdF5ELYVSqURWVhbee+89AMDIkSMxffp0dmsZqTVr1kiJzsSJE6HRaJjoENVBb8nO8ePHoVKpEBAQIC0LDAyESqWSBkKry8mTJ5GamoqZM2dWey88PBxOTk7o378/1q5dC61WW+t2li9fLnWlqVQqnRv5iFqKBQsWSD80tmzZAjMzM9y7d0/mqKi+Kisr0bVrV7zyyisAgAMHDiA6OprdVkT1oLdPSXZ2do2/NpydnaXL6XWJiopCt27dMHDgQJ3l77zzDnbu3IkjR45g6tSpWLBgAZYtW1brdt544w3k5+dLr8zMzIY1hshEdO3aFRUVFejSpQsAoG3btvj2229ljorq8ssvv8DKygoXLlwAcH+uq6eeekrmqIiMR4OTncjIyFpvIn7w+umnnwDcv7fmt4QQNS7/rdLSUmzdurXGqzp/+9vfEBQUhN69e2PBggVYunQp/vnPf9a6LaVSCXt7e50XUUtlaWmJX375BZ999hkAYNy4cZg4ceIjr46SfN5//314e3sDAKZPnw6tVos2bdrIHBWRcWnwODvh4eGYOnXqI+t07NgRp0+fxq1bt6q9l5OTAxcXlzr3s2vXLpSUlOCPf/xjnXUDAwNRUFCAW7du1WvbRAS8/PLLGDNmDDp16oSvv/4a5ubmvNHVgJSXl8PT01P6HuWN5USN1+Bkx8nJSZpB91GCgoKQn5+PEydOYMCAAQCApKQk5OfnV+uWqklUVBQmTJiAdu3a1Vk3JSUF1tbWcHBwqLsBRCTx8vJCZWUlAgICcOrUKbi4uGDXrl34/e9/L3doLdrZs2fRs2dPqZyXlweVSiVjRETGTW/37HTr1g1jxozB7NmzkZiYiMTERMyePRvjx4+XLsneuHEDPj4+OHHihM66v/76K44ePYpZs2ZV2+7+/fuxfv16nD17FhcvXsSGDRvw5ptv4oUXXoBSqdRXc4hMloWFBU6ePImNGzcCACZNmoSRI0eyW0smS5culRKdl156CVqtlokO0WPS63QRW7ZswauvvoqQkBAAwIQJE7B69Wrp/crKSmRkZKCkpERnvS+++ALt27eX1nuYpaUlPvvsM8yfPx9arRadOnXC0qVLMWfOHH02hcjkhYWFYeTIkXB3d8f3338Pc3NzTjvQjEpLS+Hk5CR9H3KaD6Kmw4lAebMykQ5OKNn8Tp06BX9/f6nMCVyJamZwgwoSkXEyNzdHfHw8tm/fDgAIDQ1FQEAAqqqqZI7MNC1atEhKdBYsWAAhBBMdoibGKzu8skNUq1u3bkGtVkvlK1euwNPTU8aITEdxcbFOUpOcnIx+/frJGBGR4eOVHSJqci4uLtBoNBg3bhyA+8NKrF+/XuaojN/x48elREepVKKkpISJDpEeMdkhokcyMzPDN998g6+//hrA/Ql9u3fvjsrKSpkjMz5CCISHh0vDb7z11lsoKyur9yTGRNQ47MZiNxZRvd29e1dnnK0LFy7giSeekDEi4/Hge+eBn3/+Gb169ZIxIiLjw24sItK7tm3bQqvVYvLkyQCALl264OOPP5Y5KsMXFxcnJTpOTk4oLS1lokPUjJjsEFGDKBQK7NixAwcPHgQAREREwMPDAxUVFTJHZniEEAgLC8OTTz4JAFi5ciVycnJgbW0tc2RELYteBxUkItMVEhKC3NxcODo6IjMzE0qlEmlpafD19ZU7NIOQm5urM2Hn+fPnpdHjiah58coOETWag4MDtFot/vznPwMAunfvjmXLlskclfxiYmKkRKdTp04oLy9nokMkIyY7RPRYFAoFoqKiEB8fDwB488034ejoiLKyMpkja35CCEyaNAljx44FAHz66ae4ePEirKysZI6MqGVjNxYRNYmhQ4eioKAA9vb2yMvLg42NDVJSUtC7d2+5Q2sWd+7cQbt27aTyxYsX0alTJxkjIqIHeGWHiJqMnZ0dtFot5s2bBwDo06cP3nzzTZmj0r+9e/dKiY6fnx8qKyuZ6BAZEI6zw3F2iPQiKSkJgYGBAO4PTFhYWIhWrVrJHFXT0mq1GDt2LA4dOgQAiIqKku5fIqKm19jzN7uxiEgvAgICUFRUBHt7e2i1Wtja2iIxMREBAQFyh9YksrKy4ObmJpWvXbsGd3d3GSMiotqwG4uI9MbW1hYajQaLFy8GAAQGBiIiIgLGfkF527ZtUqIzePBgVFVVMdEhMmDsxmI3FlGzSElJQd++faVyQUEB7OzsZIyo4TQaDYYNG4b//ve/AICtW7di2rRpMkdF1HKwG4uIDFqfPn1QUlICV1dX6YsqPj4eQ4cOlTu0esnMzISHh4dUzsrKglqtljEiIqovdmMRUbOxsbFBXl4e/vGPfwAAhg0bhlmzZhl8t9YXX3whJTqjR4+GRqNhokNkRNiNxW4sIlmcO3cO3bt3l8q5ublwcHCQMaLqqqqq4O/vj9OnTwMAoqOjMXHiRJmjImq5OOs5ERkVX19flJWVSTf2Ojo64vDhwzJH9f8uXboES0tLKdHJyclhokNkpJjsEJFslEolrl27hg8//BDA/clFp06dKnu31urVq9G5c2cAwKRJk6DVauHk5CRrTETUeOzGYjcWkUH49ddf0aVLF6l8584dtG3btlq9hIQEFBYWYvTo0U0eQ0VFBby9vXHlyhUAwHfffYcxY8Y0+X6IqHHYjUVERu2JJ55ARUUFunXrBgBwcnLC/v37pfc1Gg0iIyMxeMgQjBkzBgcPHmzS/Z8/fx5KpVJKdO7du8dEh8hEMNkhIoNhaWmJc+fO4fPPPwcATJgwAePHj8eNGzfw5PARWLr0HdgPnIZWXn0w409/Rl5eXpPsd9WqVVKSFRYWBq1WC0dHxybZNhHJj91Y7MYiMkhXr15Fx44dpbKVvRMcxy2AtUdPVBXk4NbGcEyf+gds3Lix0ft4cIP0nTt3AACxsbEIDg5+zMiJSF/YjUVEJqV9+/b461//KpWdZ3wMa4+eAAAL+3ZQPTkTmzZtwoEDB2pcf82aNejVqxfs7e1hb2+PoKAgfPfdd9L7p0+fho2NjZToFBQUMNEhMlFMdojI4GRmZmLI0GF47/0P4BAcBo/X98G8lUqnjm3PUWjVuR/+PHMW7t27V20bHTp0wIoVK/DTTz/hp59+wvDhw/HMM88gLS0NS5YsgZ+fHwAgPDwcWq3W6KauIKL6YzcWu7GIDEpubi68OnVGYUkp2k2KhLV7j1rrVhXewa2N4Zj87DPY8uWXdW67TZs2KCgogEajAXD/ya6goKAmi52I9IvdWERkEuzs7KB2dYW2ogy53/wTubFfoCLnao11LeycoBo+G1u3bMHevXsBAGWVGuQUlqOsUiPV02g0WLZsGXJzc6VEp6ioiIkOUQvBiUCJyKBYWFggPe0sTp06hU2bNuHLLVuRdWIPbNy6wNp3OGy7DdXp0rLtPhxlvyQgbOZsPHOzNX7MLIdWAGYKwN+uEAfenYWy0hJpoMKFCxdixYoVcjWPiGSg1ys77777LgYOHIhWrVrVe84bIQQiIyPh5uYGGxsbBAcHIy0tTadObm4uQkNDoVKpoFKpEBoa2mSPoBKR/BQKBfz9/fHpp5/iVnYW9uzZg5H9uqEgLgo3P5uBO9HvouRCIoSmCgqFAo4h4SgsKcfu1e9A+7+Oea0AknOA0pJiKdFxcHDAH//4RxlbRkRy0GuyU1FRgT/84Q94+eWX673OqlWr8MEHH2D16tVITk6GWq3GqFGjUFhYKNV57rnnkJqaipiYGMTExCA1NRWhoaH6aAIRyczKygq/+93vsO/rr5F18yY+/OB9dG5Vhpw9/0DWmhm4d+RzVBXdRZtRL6I4/SiKzx8DAJRdP4erH00DACisbPBj+nX4+/vj448/lrM5RCSDZrlBedOmTYiIiKjz6osQAm5uboiIiMDChQsBAOXl5XBxccHKlSvx4osvIj09Hb6+vkhMTERAQAAAIDExEUFBQTh//jy8vb3rjIc3KBMZvzNnzmDz5s341/qNKCu4B0snT1QV3AbMLNCqc38Up/0AAFANfh5thzyHUb4uuLDxdbi7u2PTpk3yBk9EjdLY87dB3bNz+fJlZGdnIyQkRFqmVCoxbNgwJCQk4MUXX8Tx48ehUqmkRAcAAgMDoVKpkJCQUGOyU15ejvLycqlcUFCg34YQkd717NkT/1i+ErvNh6Lk0ikUnf0Blbk3gYpSKdFp+9RfYO3RA6W3LmN7NsTPAgAADqVJREFU3GYUJsUhJiZG5siJqLkZVLKTnZ0NAHBxcdFZ7uLigqtXr0p1nJ2dq63r7Owsrf9by5cvx9tvv93E0RKR3ArLqiAU5rDp3B82nftDU1aEgqQ9qMq7CYWFEnnHtkBTfA9mSltYteuI7Xv2YdSoUXKHTUTNrMHJTmRkZJ2JQ3JyMvr169fooBQKhU5ZCKGz7Lfv11TnYW+88Qbmz58vlQsKCuDu7t7o+IjIMNhZW8BMAemmZHPr1nAcVvMNyGYK4OmnOLEnUUvU4GQnPDwcU6dOfWSdh+ezaQi1Wg3g/tUbV1dXafnt27elqz1qtRq3bt2qtm5OTk61K0IPKJVKKJXKRsVERIbL2tIco3xdcCT9NjTa2m8/NDdTYJSvC6wtzZsxOiIyFA1OdpycnODk5KSPWODl5QW1Wo3Dhw+jT58+AO4/0RUfH4+VK1cCAIKCgpCfn48TJ05gwIABAICkpCTk5+dj4MCBeomLiAzXrCGdcCit+g+gh2m1ArMGezVTRERkaPT66Pm1a9eQmpqKa9euQaPRIDU1FampqSgqKpLq+Pj4IDo6GsD97qmIiAgsW7YM0dHROHv2LMLCwtCqVSs899xzAIBu3bphzJgxmD17NhITE5GYmIjZs2dj/Pjx9XoSi4hMS/+ObfDOxB5Q4P4VnIeZmymgAPDOxB7o17GNLPERkfz0eoPyW2+9hc2bN0vlB1drYmNjpdmFMzIykJ+fL9V5/fXXUVpaildeeQW5ubkICAjAoUOHdCbp27JlC1599VXpqa0JEyZg9erV+mwKERmw6YGe8FHbYcOxyziUli2NoDzK1wWzBnsx0SFq4TgRKMfZITIpZZUaFJZVwc7agvfoEJkYkxhnh4jocVlbmjPJISIdnPWciIiITBqTHSIiIjJpTHaIiIjIpDHZISIiIpPGZIeIiIhMGpMdIiIiMmlMdoiIiMikMdkhIiIik8Zkh4iIiExaixxB+cEMGQUFBTJHQkRERPX14Lzd0JmuWmSyU1hYCABwd3eXORIiIiJqqMLCQqhUqnrXb5ETgWq1Wty8eRN2dnZQKBRNvv2CggK4u7sjMzPTJCcaZfuMG9tn3Ng+48b2PR4hBAoLC+Hm5gYzs/rfidMir+yYmZmhQ4cOet+Pvb29Sf4xP8D2GTe2z7ixfcaN7Wu8hlzReYA3KBMREZFJY7JDREREJs08MjIyUu4gTJG5uTmCg4NhYWGaPYVsn3Fj+4wb22fc2L7m1yJvUCYiIqKWg91YREREZNKY7BAREZFJY7JDREREJo3JDhEREZk0JjtERERk0pjs1NNnn30GLy8vWFtbw9/fHz/++OMj6+/evRu+vr5QKpXw9fVFdHS0zvtCCERGRsLNzQ02NjYIDg5GWlqaPpvwSA1p3/r16zFkyBA4OjrC0dERI0eOxIkTJ3TqhIWFQaFQ6LwCAwP13YxaNaR9mzZtqha7QqFAWVlZo7epbw2JJTg4uMb2jRs3TqpjKMfv6NGjePrpp+Hm5gaFQoG9e/fWuU58fDz8/f1hbW2NTp06Ye3atdXqGMqxa2j79uzZg1GjRqFdu3awt7dHUFAQDh48qFMnMjKy2rFTq9X6bEatGtq+uLi4Gv82z58/r1Ovru/X5tLQ9tX0uVIoFOjevbtUx5CO3/Lly9G/f3/Y2dnB2dkZEydOREZGRp3rGeL5j8lOPezYsQMRERF48803kZKSgiFDhmDs2LG4du1ajfWPHz+OKVOmIDQ0FD///DNCQ0MxefJkJCUlSXVWrVqFDz74AKtXr0ZycjLUajVGjRolTVLanBravri4OEybNg2xsbE4fvw4PDw8EBISghs3bujUGzNmDLKysqTXt99+2xzNqaah7QPuD3X+cOxZWVmwtrZ+rG3qS0Nj2bNnj067zp49C3Nzc/zhD3/QqWcIx6+4uBh+fn5YvXp1vepfvnwZTz31FIYMGYKUlBQsXrwYr776Knbv3i3VMaRj19D2HT16FKNGjcK3336LkydP4sknn8TTTz+NlJQUnXrdu3fXOXZnzpzRR/h1amj7HsjIyNCJv0uXLtJ79fl+bS4Nbd/HH3+s067MzEy0adOm2mfPUI5ffHw85syZg8TERBw+fBhVVVUICQlBcXFxresY7PlPUJ0GDBggXnrpJZ1lPj4+YtGiRTXWnzx5shgzZozOstGjR4upU6cKIYTQarVCrVaLFStWSO+XlZUJlUol1q5d28TR162h7futqqoqYWdnJzZv3iwtmzFjhnjmmWeaNM7Gamj7Nm7cKFQqVZNuU58eN5YPP/xQ2NnZiaKiImmZIR2/BwCI6OjoR9Z5/fXXhY+Pj86yF198UQQGBkplQzp2D6tP+2ri6+sr3n77bam8ZMkS4efn15ShNYn6tC82NlYAELm5ubXWqev7VS6NOX7R0dFCoVCIK1euSMsM9fgJIcTt27cFABEfH19rHUM9//HKTh0qKipw8uRJhISE6CwPCQlBQkJCjescP368Wv3Ro0dL9S9fvozs7GydOkqlEsOGDat1m/rSmPb9VklJCSorK9GmTRud5XFxcXB2dkbXrl0xe/Zs3L59u8nirq/Gtq+oqAienp7o0KEDxo8fr/PLuSn+nzWVpoglKioKU6dOha2trc5yQzh+DVXbZ++nn35CZWWlQR27pqDValFYWFjts3fhwgW4ubnBy8sLU6dOxaVLl2SKsHH69OkDV1dXjBgxArGxsTrv1fX9akyioqIwcuRIeHp66iw31OOXn58PANX+3h5mqOc/Jjt1uHPnDjQaDVxcXHSWu7i4IDs7u8Z1srOzH1n/wX8bsk19aUz7fmvRokVo3749Ro4cKS0bO3YstmzZgh9++AHvv/8+kpOTMXz4cJSXlzdp/HVpTPt8fHywadMm7Nu3D9u2bYO1tTUGDRqECxcuNHqb+vK4sZw4cQJnz57FrFmzdJYbyvFrqNo+e1VVVbhz545BHbum8P7776O4uBiTJ0+WlgUEBODf//43Dh48iPXr1yM7OxsDBw7E3bt3ZYy0flxdXbFu3Trs3r0be/bsgbe3N0aMGIGjR49Kder6fjUWWVlZ+O6776p99gz1+AkhMH/+fAwePBg9evSotZ6hnv8MZ+IKA6dQKHTKQohqyxpav6Hb1KfGxrJq1Sps27YNcXFxOve0TJkyRfp3jx490K9fP3h6euLAgQN49tlnmy7wempI+wIDA3Vuxh00aBD69u2LTz/9FJ988kmjtqlvjY0lKioKPXr0wIABA3SWG9rxa4ia/l88WP7wv39bR65j11jbtm1DZGQkvv76azg7O0vLx44dK/27Z8+eCAoKQufOnbF582bMnz9fjlDrzdvbG97e3lI5KCgImZmZeO+99zB06FBpuSkcv02bNsHBwQETJ07UWW6oxy88PBynT5/GsWPH6qxriOc/Xtmpg5OTE8zNzatlnLdv366WmT6gVqsfWf/BnfUN2aa+NKZ9D7z33ntYtmwZDh06hF69ej2yrqurKzw9PaWrI83lcdr3gJmZGfr37y/F3hTbbCqPE0tJSQm2b99e7ZdlTeQ6fg1V22fPwsICbdu2Nahj9zh27NiBmTNn4quvvtK5oloTW1tb9OzZ0+CPXW0CAwN1Yq/r+9UYCCHwxRdfIDQ0FFZWVo+sawjHb+7cudi3bx9iY2PRoUOHR9Y11PMfk506WFlZwd/fH4cPH9ZZfvjwYQwcOLDGdYKCgqrVP3TokFTfy8sLarVap05FRQXi4+Nr3aa+NKZ9APDPf/4T77zzDmJiYtCvX78693P37l1kZmbC1dX1sWNuiMa272FCCKSmpkqxN8U2m8rjxPLVV1+hvLwc06dPr3M/ch2/hqrts9evXz9YWloa1LFrrG3btiEsLAxbt27VGS6gNuXl5UhPTzf4Y1eblJQUndjr+n41BvHx8fj1118xc+bMOuvKefyEEAgPD8eePXvwww8/wMvLq851DPb8p7dbn03I9u3bhaWlpYiKihLnzp0TERERwtbWVrqDPjQ0VOdJjv/+97/C3NxcrFixQqSnp4sVK1YICwsLkZiYKNVZsWKFUKlUYs+ePeLMmTNi2rRpwtXVVRQUFBh8+1auXCmsrKzErl27RFZWlvQqLCwUQghRWFgoFixYIBISEsTly5dFbGysCAoKEu3btzeK9kVGRoqYmBhx8eJFkZKSIv70pz8JCwsLkZSUVO9tGnL7Hhg8eLCYMmVKteWGdPwKCwtFSkqKSElJEQDEBx98IFJSUsTVq1eFEEIsWrRIhIaGSvUvXbokWrVqJf7yl7+Ic+fOiaioKGFpaSl27dol1TGkY9fQ9m3dulVYWFiIf/3rXzqfvby8PKnOggULRFxcnLh06ZJITEwU48ePF3Z2dkbRvg8//FBER0eLX375RZw9e1YsWrRIABC7d++W6tTn+9VQ2/fA9OnTRUBAQI3bNKTj9/LLLwuVSiXi4uJ0/t5KSkqkOsZy/mOyU0//+te/hKenp7CyshJ9+/bVefRu2LBhYsaMGTr1d+7cKby9vYWlpaXw8fHR+bAKcf/xuyVLlgi1Wi2USqUYOnSoOHPmTHM0pUYNaZ+np6cAUO21ZMkSIYQQJSUlIiQkRLRr105YWloKDw8PMWPGDHHt2rVmbtX/a0j7IiIihIeHh7CyshLt2rUTISEhIiEhoUHbbG4N/fvMyMgQAMShQ4eqbcuQjt+DR5F/+3rQnhkzZohhw4bprBMXFyf69OkjrKysRMeOHcWaNWuqbddQjl1D2zds2LBH1hdCiClTpghXV1dhaWkp3NzcxLPPPivS0tKat2H/09D2rVy5UnTu3FlYW1sLR0dHMXjwYHHgwIFq263r+7W5NObvMy8vT9jY2Ih169bVuE1DOn41tQ2A2Lhxo1THWM5/iv81iIiIiMgk8Z4dIiIiMmlMdoiIiMikMdkhIiIik8Zkh4iIiEwakx0iIiIyaUx2iIiIyKQx2SEiIiKTxmSHiIiITBqTHSIiIjJpTHaIiIjIpDHZISIiIpP2f3Z2MkhX7/qbAAAAAElFTkSuQmCC",
"text/plain": [
"Figure(PyObject )"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"#This function represent a Graph\n",
"function plot_graph(G,plot_edges = false, plot_vertices = false)\n",
" posx = zeros(nv(G))\n",
" posy = zeros(nv(G))\n",
" for v in 1:nv(G)\n",
" posx[v] = get_prop(G,v,:x)\n",
" posy[v] = get_prop(G,v,:y)\n",
" if plot_vertices\n",
" text(posx[v],posy[v],v)\n",
" end\n",
" end\n",
" scatter(posx, posy, s=50) # plotting edges\n",
" \n",
" for e in edges(G) #edge s -> t\n",
" x_s,y_s = posx[src(e)], posy[src(e)]\n",
" x_t,y_t = posx[dst(e)], posy[dst(e)]\n",
" arrow(x_s, y_s, x_t - x_s, y_t - y_s,\n",
" head_width=0.05, length_includes_head=true)\n",
" if plot_edges\n",
" δ_x, δ_y = 0.05.*(y_s-y_t , x_t - x_s)\n",
" text(0.5*x_s+0.5*x_t + δ_x, 0.5*y_s+0.5*y_t + δ_y, get_prop(G,e,:idx)) # edges number\n",
" end\n",
" end\n",
"end\n",
"plot_graph(G,true,true)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"matrix_rep (generic function with 1 method)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Matrix representation of conservation law :: Nx = io\n",
"function matrix_rep(G)\n",
" io = [get_prop(G,i,:io) for i=1:nv(G)]\n",
"\n",
" N = zeros(Int,(nv(G),ne(G)))\n",
" for e in edges(G)\n",
" j = get_prop(G,e,:idx)\n",
" N[src(e),j] = -1 \n",
" N[dst(e),j] = 1\n",
" end\n",
" \n",
" return N,io\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Computing the price of anarchy\n",
"\n",
"Define $C(x)=\\sum_{e\\in E}x_e \\ell_e(x_e)$ the price associated with the \n",
"admissible flux $x$. We denote by $x^{UE}$ the user equilibrium, and $x^{SO}$\n",
"the social optimum. The price of anarchy is by definition \n",
"$$PoA = \\frac{C(x^{UE})}{C(x^{SO})} \\geq 1 .$$"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"social_cost (generic function with 1 method)"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"### Computing social cost associated to an intensity vector\n",
"function social_cost(G,x)\n",
" cost = 0\n",
" for e in edges(G)\n",
" idx = get_prop(G,e,:idx) # getting index of edge \n",
" a = get_prop(G,e,:a) # get coefficient a of the cost function\n",
" b = get_prop(G,e,:b) # get coefficient b of the cost function\n",
" cost += a*x[idx]^2 + b*x[idx] \n",
" end\n",
" return cost\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"******************************************************************************\n",
"This program contains Ipopt, a library for large-scale nonlinear optimization.\n",
" Ipopt is released as open source code under the Eclipse Public License (EPL).\n",
" For more information visit http://projects.coin-or.org/Ipopt\n",
"******************************************************************************\n",
"\n",
"Edge 1 => 2 has intensity 0.5\n",
"Edge 1 => 3 has intensity 0.5\n",
"Edge 2 => 4 has intensity 0.5\n",
"Edge 3 => 4 has intensity 0.5\n"
]
},
{
"data": {
"text/plain": [
"1.5"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"### Computing Social Optimum for graph G\n",
"function SO(G; display = false)\n",
" N,io = matrix_rep(G)\n",
" # Setting up the optimization model\n",
" m_so = Model(optimizer_with_attributes(Ipopt.Optimizer, \"print_level\" => 0))\n",
" \n",
" @variable(m_so,x_so[1:ne(G)] >= 0) # (non-negative) edge-intensity\n",
" @variable(m_so,cost_so[1:ne(G)]) # cost by edge\n",
" @constraint(m_so,N*x_so .== io) # Conservation law\n",
" \n",
"\n",
" # Defining social optimum cost \n",
" for e in edges(G)\n",
" idx = get_prop(G,e,:idx) # getting index of edge \n",
" a = get_prop(G,e,:a) # get coefficient a of the cost function\n",
" b = get_prop(G,e,:b) # get coefficient b of the cost function\n",
" @NLconstraint(m_so,cost_so[idx] == a*x_so[idx]^2 + b*x_so[idx]) \n",
" end\n",
" @NLobjective(m_so,Min,sum(cost_so[idx] for idx=1:ne(G)))\n",
"\n",
" # Solving the optimisation problem\n",
" optimize!(m_so)\n",
" \n",
" x_opt = zeros(ne(G))\n",
" for e in edges(G)\n",
" x_opt[get_prop(G,e,:idx)] = round(value(x_so[get_prop(G,e,:idx)]),digits=3)\n",
" if display\n",
" # Printing the computed intensity on all edges\n",
" print(e)\n",
" print(\" has intensity \")\n",
" println(x_opt[get_prop(G,e,:idx)])\n",
" end\n",
" end\n",
" \n",
" return x_opt # return the optimal intensity\n",
"end\n",
"\n",
"x_so = SO(G, display=true)\n",
"\n",
"round(social_cost(G,x_so),digits=3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 6\n",
"\n",
"From the function computing the social optimum, define a function computing the user equilibrium. \n",
"\n",
"Deduce the price of anarchy of this network. Comment the solution."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Edge 1 => 2 has intensity 0.5\n",
"Edge 1 => 3 has intensity 0.5\n",
"Edge 2 => 4 has intensity 0.5\n",
"Edge 3 => 4 has intensity 0.5\n"
]
},
{
"data": {
"text/plain": [
"1.5"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"### Computing Social Optimum for graph G\n",
"function UE(G; display = false)\n",
" N,io = matrix_rep(G)\n",
" # Setting up the optimization model\n",
" m_ue = Model(optimizer_with_attributes(Ipopt.Optimizer, \"print_level\" => 0))\n",
" \n",
" @variable(m_ue,x_ue[1:ne(G)] >= 0) # (non-negative) edge-intensity\n",
" @variable(m_ue,cost_ue[1:ne(G)]) # cost by edge\n",
" @constraint(m_ue,N*x_ue .== io) # Conservation law\n",
" \n",
"\n",
" # Defining social optimum cost \n",
" for e in edges(G)\n",
" idx = get_prop(G,e,:idx) # getting index of edge \n",
" a = get_prop(G,e,:a) # get coefficient a of the cost function\n",
" b = get_prop(G,e,:b) # get coefficient b of the cost function\n",
" # \n",
" # A COMPLETER \n",
" #@NLconstraint(m_ue,cost_ue[idx] == ????? )\n",
" #\n",
" #\n",
" end\n",
" @NLobjective(m_ue,Min,sum(cost_ue[idx] for idx=1:ne(G)))\n",
"\n",
" # Solving the optimisation problem\n",
" optimize!(m_ue)\n",
" \n",
" x_opt = zeros(ne(G))\n",
" for e in edges(G)\n",
" x_opt[get_prop(G,e,:idx)] = round(value(x_ue[get_prop(G,e,:idx)]),digits=3)\n",
" if display\n",
" # Printing the computed intensity on all edges\n",
" print(e)\n",
" print(\" has intensity \")\n",
" println(x_opt[get_prop(G,e,:idx)])\n",
" end\n",
" end\n",
" \n",
" return x_opt # return the optimal intensity\n",
"end\n",
"\n",
"x_ue = UE(G, display=true)\n",
"\n",
"round(social_cost(G,x_ue),digits=3)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.0"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function PoA(G)\n",
" x_so = SO(G)\n",
" x_ue = UE(G)\n",
" # \n",
" #\n",
" # A COMPLETER \n",
" # return ????\n",
"end\n",
"PoA(G)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 7\n",
"\n",
"Add one edge form node $2$ to $3$ with null cost. Compute the price of anarchy.\n",
"\n",
"Add one edge form node $3$ to $2$ with null cost. Compute the price of anarchy.\n",
"\n",
"Add both edges. Compute the price of anarchy.\n",
"\n",
"Look at obtained intensities, and comment the result."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.3333333333333333"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# \n",
"# A COMPLETER \n",
"# "
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.0"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# \n",
"# A COMPLETER \n",
"# "
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.3333333333333333"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# \n",
"# A COMPLETER \n",
"# "
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"6-element Array{Float64,1}:\n",
" 1.0\n",
" 0.0\n",
" 0.0\n",
" 1.0\n",
" 100000.5\n",
" 99999.5"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# \n",
"# A COMPLETER \n",
"# "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 8\n",
"\n",
"In the previous graph modify the constant cost functions equal to $1$\n",
"into cost of the form $1+0.5x$. What do you observe on optimal solutions\n",
"$x^{UE}$ and $x^{S0}$ ? Comment."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{4, 6} directed Int64 metagraph with Float64 weights defined by :weight (default weight 1.0)"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# \n",
"# A COMPLETER \n",
"# "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 9\n",
"\n",
"Choose random affine latency on each edge, the coefficients being chosen uniformly on\n",
"$[0,1]$ (using rand()). Plot an histogram of the price of anarchy over\n",
"$100$ realization. What result from the course do we observe ?"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAhoAAAGgCAYAAADsAM6oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAHeJJREFUeJzt3X+Q1PV9+PHXeRyLGiDxED10UcEJIGibQmqI+CsxWqrU/NFkTMVSk7RpvSTg2TZcrFES4TRpLJ0mwUoZaocqTmKwjmOMOg2oMTRCobV6QQ3qYdXYJcqCwMqPz/cPx/vmIhg+y76X2/PxmNmZfDafz2dfvu+ze8+5u2WbsizLAgAggcMO9QAAwMAlNACAZIQGAJCM0AAAkhEaAEAyQgMASEZoAADJCA0AIBmhAQAkIzQAgGSEBgCQzKB6P+DevXvjxRdfjKFDh0ZTU1O9Hx4AqEKWZbF169YYNWpUHHbYgf+cou6h8eKLL0axWKz3wwIANbBp06Y4/vjjD3j/uofG0KFDI+LNQYcNG1bvhwcAqlAul6NYLPZ+Hz9QdQ+Nt35dMmzYMKEBAA0m7589+GNQACAZoQEAJCM0AIBkhAYAkIzQAACSERoAQDJCAwBIRmgAAMkIDQAgGaEBACQjNACAZHKFxu7du+Nv/uZv4qSTTorDDz88xowZE1/96ldj7969qeYDABpYrg9Vu/HGG+Pmm2+OW2+9NSZOnBhr1qyJyy+/PIYPHx6zZ89ONSMA0KByhcZPfvKTuPjii+PCCy+MiIgTTzwxbr/99lizZs1+j6lUKlGpVHq3y+VylaMOTD09PVEqlZKce8SIETF69Ogk504l1Xo04lpQP647SCdXaEybNi1uvvnmeOqpp+L9739//Nd//Vc88sgjsXDhwv0e09XVFfPmzTvoQQeinp6eGDd+QuzcsT3J+YccfkRs+Fl3w7zQpVyPRlsL6sd1B2nlCo0vfelLsWXLlhg/fnw0NzfHnj17Yv78+fGpT31qv8d0dnZGR0dH73a5XI5isVj9xANIqVSKnTu2R+tFV0VLa23XZNfmTbH5nm9GqVRqmBe5VOvRiGtB/bjuIK1coXHHHXfEsmXL4rbbbouJEyfG+vXrY86cOTFq1KiYNWvWPo8pFApRKBRqMuxA1dJajMKxJx/qMfoN68Gh4LqDNHKFxl/91V/F3Llz45JLLomIiFNPPTWef/756Orq2m9oAADvXrne3rp9+/Y47LC+hzQ3N3t7KwCwT7l+ojFjxoyYP39+jB49OiZOnBjr1q2Lm266KT796U+nmg8AaGC5QuMf/uEf4pprrokrrrgiXnnllRg1alR87nOfi6985Sup5gMAGliu0Bg6dGgsXLjwHd/OCgDwFp91AgAkIzQAgGSEBgCQjNAAAJIRGgBAMkIDAEhGaAAAyQgNACAZoQEAJCM0AIBkhAYAkIzQAACSERoAQDJCAwBIRmgAAMkIDQAgGaEBACQjNACAZIQGAJCM0AAAkhEaAEAyQgMASEZoAADJCA0AIBmhAQAkIzQAgGSEBgCQjNAAAJIRGgBAMkIDAEhGaAAAyQgNACAZoQEAJJMrNE488cRoamp62629vT3VfABAAxuUZ+fHHnss9uzZ07v9P//zP/Gxj30sPvGJT9R8MACg8eUKjaOPPrrP9g033BBjx46Ns88+e7/HVCqVqFQqvdvlcjnniP1DT09PlEqlmp6zu7u7pueD/iDFcyUiYsSIETF69Oian5fGl+qai3Dd1UKu0PhVb7zxRixbtiw6Ojqiqalpv/t1dXXFvHnzqn2YfqGnpyfGjZ8QO3dsP9SjQL+W8rky5PAjYsPPur3o00fq12fX3cGrOjTuuuuueO211+JP/uRP3nG/zs7O6Ojo6N0ul8tRLBarfdhDolQqxc4d26P1oquipbV2s+/YuCa2PLysZueDQy3Vc2XX5k2x+Z5vRqlU8oJPH6muuQjXXa1UHRpLliyJ6dOnx6hRo95xv0KhEIVCodqH6VdaWotROPbkmp1v1+ZNNTsX9Ce1fq7Ab+Ka67+qCo3nn38+Hnzwwfj+979f63kAgAGkqn9HY+nSpTFy5Mi48MILaz0PADCA5A6NvXv3xtKlS2PWrFkxaFDVv3kBAN4FcofGgw8+GD09PfHpT386xTwAwACS+0cS559/fmRZlmIWAGCA8VknAEAyQgMASEZoAADJCA0AIBmhAQAkIzQAgGSEBgCQjNAAAJIRGgBAMkIDAEhGaAAAyQgNACAZoQEAJCM0AIBkhAYAkIzQAACSERoAQDJCAwBIRmgAAMkIDQAgGaEBACQjNACAZIQGAJCM0AAAkhEaAEAyQgMASEZoAADJCA0AIBmhAQAkIzQAgGSEBgCQjNAAAJIRGgBAMrlD43//939j5syZ0draGkcccUT89m//dqxduzbFbABAgxuUZ+dXX301zjjjjDj33HPjBz/4QYwcOTJ+/vOfx3vf+95U8wEADSxXaNx4441RLBZj6dKlvfedeOKJ73hMpVKJSqXSu10ul/NNmENPT0+USqWan7e7u7vm5wTyS/Fc9PyunxSv0b5+/V+u0Lj77rvjggsuiE984hOxatWqOO644+KKK66IP/3TP93vMV1dXTFv3ryDHvQ36enpiXHjJ8TOHduTPxZQX3u2vRrR1BQzZ8481KNQJa/R7165QmPjxo2xaNGi6OjoiC9/+cvx05/+NL74xS9GoVCIP/7jP97nMZ2dndHR0dG7XS6Xo1gsHtzU+1AqlWLnju3RetFV0dJa2/Pv2Lgmtjy8rKbnBA7c3sq2iCzz/G5gqV6jff36v1yhsXfv3pgyZUosWLAgIiI+8IEPxBNPPBGLFi3ab2gUCoUoFAoHP+kBamktRuHYk2t6zl2bN9X0fEB1PL8bX62/hr5+/V+ud520tbXFKaec0ue+CRMmRE9PT02HAgAGhlyhccYZZ8SGDRv63PfUU0/FCSecUNOhAICBIVdoXHnllbF69epYsGBBPPPMM3HbbbfFLbfcEu3t7anmAwAaWK7Q+OAHPxgrVqyI22+/PSZNmhRf+9rXYuHChXHppZemmg8AaGC5/hg0IuKiiy6Kiy66KMUsAMAA47NOAIBkhAYAkIzQAACSERoAQDJCAwBIRmgAAMkIDQAgGaEBACQjNACAZIQGAJCM0AAAkhEaAEAyQgMASEZoAADJCA0AIBmhAQAkIzQAgGSEBgCQjNAAAJIRGgBAMkIDAEhGaAAAyQgNACAZoQEAJCM0AIBkhAYAkIzQAACSERoAQDJCAwBIRmgAAMkIDQAgGaEBACQjNACAZHKFxnXXXRdNTU19bscee2yq2QCABjco7wETJ06MBx98sHe7ubm5pgMBAANH7tAYNGhQrp9iVCqVqFQqvdvlcjnvQ0K/1tPTE6VSKcm5R4wYEaNHj675eVPN3N3dXfNzwkCV6nmY6nWjWrlD4+mnn45Ro0ZFoVCI008/PRYsWBBjxozZ7/5dXV0xb968gxoS+quenp4YN35C7NyxPcn5hxx+RGz4WXdNXzRSzwz8ZimfhyleNw5GrtA4/fTT41/+5V/i/e9/f/ziF7+I66+/Pj784Q/HE088Ea2trfs8prOzMzo6Onq3y+VyFIvFg5sa+olSqRQ7d2yP1ouuipbW2l7XuzZvis33fDNKpVJNXzBSzrxj45rY8vCymp4TBqJUz8NUrxsHI1doTJ8+vfd/n3rqqTF16tQYO3Zs3HrrrX1i4lcVCoUoFAoHNyX0cy2txSgce/KhHiOXFDPv2ryppueDga4RXzvyOqi3tx555JFx6qmnxtNPP12reQCAAeSgQqNSqUR3d3e0tbXVah4AYADJFRp/+Zd/GatWrYpnn302/uM//iP+8A//MMrlcsyaNSvVfABAA8v1NxovvPBCfOpTn4pSqRRHH310fOhDH4rVq1fHCSeckGo+AKCB5QqN5cuXp5oDABiAfNYJAJCM0AAAkhEaAEAyQgMASEZoAADJCA0AIBmhAQAkIzQAgGSEBgCQjNAAAJIRGgBAMkIDAEhGaAAAyQgNACAZoQEAJCM0AIBkhAYAkIzQAACSERoAQDJCAwBIRmgAAMkIDQAgGaEBACQjNACAZIQGAJCM0AAAkhEaAEAyQgMASEZoAADJCA0AIBmhAQAkIzQAgGSEBgCQzEGFRldXVzQ1NcWcOXNqNQ8AMIBUHRqPPfZY3HLLLXHaaafVch4AYAAZVM1B27Zti0svvTQWL14c119//TvuW6lUolKp9G6Xy+VqHpJ+pKenJ0qlUs3P293dXfNzDgS1XhfrDNRTVaHR3t4eF154YZx33nm/MTS6urpi3rx5VQ1H/9PT0xPjxk+InTu2H+pRBrw9216NaGqKmTNnHupRAKqWOzSWL18e//mf/xmPPfbYAe3f2dkZHR0dvdvlcjmKxWLeh6WfKJVKsXPH9mi96Kpoaa3t13HHxjWx5eFlNT1nI9tb2RaRZTVfa+sM1FOu0Ni0aVPMnj077r///hgyZMgBHVMoFKJQKFQ1HP1XS2sxCseeXNNz7tq8qabnGyhqvdbWGainXKGxdu3aeOWVV2Ly5Mm99+3Zsyceeuih+Na3vhWVSiWam5trPiQA0JhyhcZHP/rRePzxx/vcd/nll8f48ePjS1/6ksgAAPrIFRpDhw6NSZMm9bnvyCOPjNbW1rfdDwDgXwYFAJKp6u2tv2rlypU1GAMAGIj8RAMASEZoAADJCA0AIBmhAQAkIzQAgGSEBgCQjNAAAJIRGgBAMkIDAEhGaAAAyQgNACAZoQEAJCM0AIBkhAYAkIzQAACSERoAQDJCAwBIRmgAAMkIDQAgGaEBACQjNACAZIQGAJCM0AAAkhEaAEAyQgMASEZoAADJCA0AIBmhAQAkIzQAgGSEBgCQjNAAAJIRGgBAMkIDAEgmV2gsWrQoTjvttBg2bFgMGzYspk6dGj/4wQ9SzQYANLhcoXH88cfHDTfcEGvWrIk1a9bERz7ykbj44ovjiSeeSDUfANDABuXZecaMGX2258+fH4sWLYrVq1fHxIkT93lMpVKJSqXSu10ul6sYk2p1d3f36/MB1enp6YlSqVTz844YMSJGjx5d8/Py7pUrNH7Vnj174rvf/W68/vrrMXXq1P3u19XVFfPmzav2YajSnm2vRjQ1xcyZMw/1KECN9fT0xLjxE2Lnju01P/eQw4+IDT/rFhvUTO7QePzxx2Pq1Kmxc+fOeM973hMrVqyIU045Zb/7d3Z2RkdHR+92uVyOYrFY3bQcsL2VbRFZFq0XXRUtrbVb7x0b18SWh5fV7HxAfqVSKXbu2F7z5/euzZti8z3fjFKpJDSomdyhMW7cuFi/fn289tprceedd8asWbNi1apV+42NQqEQhULhoAelOi2txSgce3LNzrdr86aanQs4OLV+fkMKuUNj8ODBcfLJb17YU6ZMicceeyz+/u//Pv7xH/+x5sMBAI3toP8djSzL+vyxJwDAW3L9ROPLX/5yTJ8+PYrFYmzdujWWL18eK1eujPvuuy/VfABAA8sVGr/4xS/isssui5deeimGDx8ep512Wtx3333xsY99LNV8AEADyxUaS5YsSTUHADAA+awTACAZoQEAJCM0AIBkhAYAkIzQAACSERoAQDJCAwBIRmgAAMkIDQAgGaEBACQjNACAZIQGAJCM0AAAkhEaAEAyQgMASEZoAADJCA0AIBmhAQAkIzQAgGSEBgCQjNAAAJIRGgBAMkIDAEhGaAAAyQgNACAZoQEAJCM0AIBkhAYAkIzQAACSERoAQDJCAwBIRmgAAMkIDQAgmVyh0dXVFR/84Adj6NChMXLkyPj4xz8eGzZsSDUbANDgcoXGqlWror29PVavXh0PPPBA7N69O84///x4/fXXU80HADSwQXl2vu+++/psL126NEaOHBlr166Ns846a5/HVCqVqFQqvdvlcrmKMeHgdXd3N8Q5GVga8bprxJlTsh4HJ1do/LotW7ZERMRRRx213326urpi3rx5B/MwcFD2bHs1oqkpZs6ceahH4V2kEa+7Rpw5JetRG1WHRpZl0dHREdOmTYtJkybtd7/Ozs7o6Ojo3S6Xy1EsFqt9WMhtb2VbRJZF60VXRUtrba+9HRvXxJaHl9X0nAwMjXjdNeLMKVmP2qg6ND7/+c/Hf//3f8cjjzzyjvsVCoUoFArVPgzUTEtrMQrHnlzTc+7avKmm52PgacTrrhFnTsl6HJyqQuMLX/hC3H333fHQQw/F8ccfX+uZAIABIldoZFkWX/jCF2LFihWxcuXKOOmkk1LNBQAMALlCo729PW677bb4t3/7txg6dGi8/PLLERExfPjwOPzww5MMCAA0rlz/jsaiRYtiy5Ytcc4550RbW1vv7Y477kg1HwDQwHL/6gQA4ED5rBMAIBmhAQAkIzQAgGSEBgCQjNAAAJIRGgBAMkIDAEhGaAAAyQgNACAZoQEAJCM0AIBkhAYAkIzQAACSERoAQDJCAwBIRmgAAMkIDQAgGaEBACQjNACAZIQGAJCM0AAAkhEaAEAyQgMASEZoAADJCA0AIBmhAQAkIzQAgGSEBgCQjNAAAJIRGgBAMkIDAEhGaAAAyQgNACCZ3KHx0EMPxYwZM2LUqFHR1NQUd911V4q5AIABIHdovP766/Fbv/Vb8a1vfSvFPADAADIo7wHTp0+P6dOnH/D+lUolKpVK73a5XM77kABAg0r+NxpdXV0xfPjw3luxWEz9kABAP5E8NDo7O2PLli29t02bNqV+SACgn8j9q5O8CoVCFAqF1A8DAPRD3t4KACQjNACAZHL/6mTbtm3xzDPP9G4/++yzsX79+jjqqKNi9OjRNR0OAGhsuUNjzZo1ce655/Zud3R0RETErFmz4p//+Z9rNhgA0Phyh8Y555wTWZalmAUAGGD8jQYAkIzQAACSERoAQDJCAwBIRmgAAMkIDQAgGaEBACQjNACAZIQGAJCM0AAAkhEaAEAyQgMASEZoAADJCA0AIBmhAQAkIzQAgGSEBgCQjNAAAJIRGgBAMkIDAEhGaAAAyQgNACAZoQEAJCM0AIBkhAYAkIzQAACSERoAQDJCAwBIRmgAAMkIDQAgGaEBACQjNACAZIQGAJBMVaHxne98J0466aQYMmRITJ48OR5++OFazwUADAC5Q+OOO+6IOXPmxNVXXx3r1q2LM888M6ZPnx49PT0p5gMAGtigvAfcdNNN8ZnPfCY++9nPRkTEwoUL44c//GEsWrQourq63rZ/pVKJSqXSu71ly5aIiCiXy9XOvE/btm178/Fefib2vrGzpufetXlTknOnOm/Kc5u58c9t5vqc28z1OXcjzpzy3Lt++UJEvPk9sdbfZ986X5Zl+Q7McqhUKllzc3P2/e9/v8/9X/ziF7Ozzjprn8dce+21WUS4ubm5ubm5DYDbpk2b8qRDlusnGqVSKfbs2RPHHHNMn/uPOeaYePnll/d5TGdnZ3R0dPRu7927N375y19Ga2trNDU15Xn4fSqXy1EsFmPTpk0xbNiwgz4f+2et68da14d1rh9rXT+p1jrLsti6dWuMGjUq13G5f3USEW8LhCzL9hsNhUIhCoVCn/ve+973VvOw72jYsGEu3jqx1vVjrevDOtePta6fFGs9fPjw3Mfk+mPQESNGRHNz89t+evHKK6+87accAAC5QmPw4MExefLkeOCBB/rc/8ADD8SHP/zhmg4GADS+5uuuu+66PAcMGzYsrrnmmjjuuONiyJAhsWDBgvjRj34US5cuTfIrkQPR3Nwc55xzTgwaVNVvgsjBWtePta4P61w/1rp++tNaN2W536fy5j/Y9fWvfz1eeumlmDRpUvzd3/1dnHXWWSnmAwAaWFWhAQBwIHzWCQCQjNAAAJIRGgBAMkIDAEimIUIj78fS33nnnXHKKadEoVCIU045JVasWFGnSRtfnrVevHhxnHnmmfG+970v3ve+98V5550XP/3pT+s4bWPLe12/Zfny5dHU1BQf//jHE084MORd59deey3a29ujra0thgwZEhMmTIh77723TtM2trxrvXDhwhg3blwcfvjhUSwW48orr4ydO2v74WUDzUMPPRQzZsyIUaNGRVNTU9x1112/8ZhVq1bF5MmTY8iQITFmzJi4+eab6zDpr8j1ySiHwPLly7OWlpZs8eLF2ZNPPpnNnj07O/LII7Pnn39+n/s/+uijWXNzc7ZgwYKsu7s7W7BgQTZo0KBs9erVdZ688eRd6z/6oz/Kvv3tb2fr1q3Luru7s8svvzwbPnx49sILL9R58saTd63f8txzz2XHHXdcduaZZ2YXX3xxnaZtXHnXuVKpZFOmTMl+//d/P3vkkUey5557Lnv44Yez9evX13nyxpN3rZctW5YVCoXsX//1X7Nnn302++EPf5i1tbVlc+bMqfPkjeXee+/Nrr766uzOO+/MIiJbsWLFO+6/cePG7Igjjshmz56dPfnkk9nixYuzlpaW7Hvf+16dJs6yfh8av/u7v5v9+Z//eZ/7xo8fn82dO3ef+3/yk5/Mfu/3fq/PfRdccEF2ySWXJJtxoMi71r9u9+7d2dChQ7Nbb701xXgDSjVrvXv37uyMM87I/umf/imbNWuW0DgAedd50aJF2ZgxY7I33nijHuMNKHnXur29PfvIRz7S576Ojo5s2rRpyWYcaA4kNP76r/86Gz9+fJ/7Pve5z2Uf+tCHUo7WR7/+1ckbb7wRa9eujfPPP7/P/eeff348+uij+zzmJz/5ydv2v+CCC/a7P2+qZq1/3fbt22PXrl1x1FFHpRhxwKh2rb/61a/G0UcfHZ/5zGdSjzggVLPOd999d0ydOjXa29vjmGOOiUmTJsWCBQtiz5499Ri5YVWz1tOmTYu1a9f2/rp148aNce+998aFF16YfN53k/19T1yzZk3s2rWrLjMc+n+b9B1U87H0L7/8cq79eVM1a/3r5s6dG8cdd1ycd955KUYcMKpZ6x//+MexZMmSWL9+fT1GHBCqWeeNGzfGv//7v8ell14a9957bzz99NPR3t4eu3fvjq985Sv1GLshVbPWl1xySfzf//1fTJs2LbIsi927d8df/MVfxNy5c+sx8rvG/r4n7t69O0qlUrS1tSWfoV+HxlvyfCx9Nfvz/1W7dl//+tfj9ttvj5UrV8aQIUNSjTegHOhab926NWbOnBmLFy+OESNG1Gu8ASPPNb13794YOXJk3HLLLdHc3ByTJ0+OF198Mb7xjW8IjQOQZ61XrlwZ8+fPj+985ztx+umnxzPPPBOzZ8+Otra2uOaaa+ox7rvGvr4u+7o/lX4dGtV8LP2xxx7rY+yrUM1av+Vv//ZvY8GCBfHggw/GaaedlnLMASHvWv/85z+P5557LmbMmNF73969eyMiYtCgQbFhw4YYO3Zs2qEbUDXXdFtbW7S0tERzc3PvfRMmTIiXX3453njjjRg8eHDSmRtVNWt9zTXXxGWXXRaf/exnIyLi1FNPjddffz3+7M/+LK6++uo47LB+/Zv9hrG/74mDBg2K1tbWuszQr7+S1Xws/dSpU9+2//333+9j7H+DatY6IuIb3/hGfO1rX4v77rsvpkyZknrMASHvWo8fPz4ef/zxWL9+fe/tD/7gD+Lcc8+N9evXR7FYrNfoDaWaa/qMM86IZ555pjfkIiKeeuqpaGtrExnvoJq13r59+9tiorm5ObI336SQbNZ3m/19T5wyZUq0tLTUZ4i6/dlpld56y9SSJUuyJ598MpszZ0525JFHZs8991yWZVl22WWX9fmr5h//+MdZc3NzdsMNN2Td3d3ZDTfc4O2tByjvWt94443Z4MGDs+9973vZSy+91HvbunXrofpPaBh51/rXedfJgcm7zj09Pdl73vOe7POf/3y2YcOG7J577slGjhyZXX/99YfqP6Fh5F3ra6+9Nhs6dGh2++23Zxs3bszuv//+bOzYsdknP/nJQ/Wf0BC2bt2arVu3Llu3bl0WEdlNN92UrVu3rvdtxHPnzs0uu+yy3v3fenvrlVdemT355JPZkiVLvL11X7797W9nJ5xwQjZ48ODsd37nd7JVq1b1/n9nn312NmvWrD77f/e7383GjRuXtbS0ZOPHj8/uvPPOOk/cuPKs9QknnJBFxNtu1157bf0Hb0B5r+tfJTQOXN51fvTRR7PTTz89KxQK2ZgxY7L58+dnu3fvrvPUjSnPWu/atSu77rrrsrFjx2ZDhgzJisVidsUVV2SvvvrqIZi8cfzoRz/a5+vuW2s7a9as7Oyzz+5zzMqVK7MPfOAD2eDBg7MTTzwxW7RoUV1n9jHxAEAy/fpvNACAxiY0AIBkhAYAkIzQAACSERoAQDJCAwBIRmgAAMkIDQAgGaEBACQjNACAZIQGAJDM/wNLrUR9hCpESQAAAABJRU5ErkJggg==",
"text/plain": [
"Figure(PyObject )"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"┌ Warning: `getindex(o::PyObject, s::Symbol)` is deprecated in favor of dot overloading (`getproperty`) so elements should now be accessed as e.g. `o.s` instead of `o[:s]`.\n",
"│ caller = top-level scope at In[15]:3\n",
"└ @ Core In[15]:3\n"
]
}
],
"source": [
"# Plotting an histogram\n",
"y = rand(100)\n",
"PyPlot.plt[:hist](y, bins=20, edgecolor=\"k\");"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGgCAYAAACXJAxkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XtwVOX9x/HPmssGMKQkkVxkoUERhACFxCo3BcVgUJTaFlChcUoZtVykgSqBdkRbCGK9DdR4YxAFGqYCSrWiQUkQEYVIkFsRFUnUxEwoZpMAC4Tz+6M/dgwEyYZz8rDJ+zVzZjjnPPvs95uFzYfds8+6LMuyBAAAYMBFpgsAAAAtF0EEAAAYQxABAADGEEQAAIAxBBEAAGAMQQQAABhDEAEAAMYQRAAAgDEEEQAAYAxBBAAAGEMQAQAAxoSaLuB0J0+e1LfffqvIyEi5XC7T5QAAgAawLEtVVVVKTEzURRc1/HWOCy6IfPvtt/J4PKbLAAAAjVBSUqIOHTo0ePwFF0QiIyMl/a+Rtm3bGq4GAAA0hNfrlcfj8f8eb6gLLoicejumbdu2BBEAAIJMoJdVBHSxak5Ojnr16uUPCf369dNbb73lPz948GC5XK4625gxYwIqCAAAtBwBvSLSoUMHzZs3T5dffrkkacmSJbrtttu0bds29ejRQ5I0YcIEPfLII/7btGrVysZyAQBAcxJQEBkxYkSd/Tlz5ignJ0ebN2/2B5HWrVsrPj7evgoBAECz1eh1RGpra5Wbm6uamhr169fPf3zZsmWKjY1Vjx49NH36dFVVVf3oPD6fT16vt84GAABahoAvVt2xY4f69euno0eP6uKLL9bq1avVvXt3SdJdd92lpKQkxcfHa+fOncrKytL27duVl5d31vmys7P18MMPN74DAAAQtFyWZVmB3ODYsWMqLi7W999/r5UrV+rFF19UQUGBP4z8UGFhoVJTU1VYWKi+ffvWO5/P55PP5/Pvn/r4T2VlJZ+aAQAgSHi9XkVFRQX8+zvgIHK6oUOH6rLLLtNzzz13xjnLsuR2u/XKK69o9OjRDZqvsY0AAABzGvv7+7y/a8ayrDqvaPzQrl27dPz4cSUkJJzv3QAAgGYooGtEZs6cqfT0dHk8HlVVVSk3N1f5+flau3atvvjiCy1btkzDhw9XbGysdu/erWnTpqlPnz4aMGCAU/UDAIAgFlAQ+e677zRu3DiVlpYqKipKvXr10tq1a3XjjTeqpKRE7777rp5++mlVV1fL4/Ho5ptv1kMPPaSQkBCn6gcAAEHsvK8RsRvXiAAAEHyMXSMCAADQWAQRAABgzAX37btOKy4uVkVFhSNzx8bGqmPHjo7MDQBAc9SigkhxcbG6drtSR48cdmT+iFattfc/ewgjAAA0UIsKIhUVFTp65LBibpmmsBiPrXMfP1iig288roqKCoIIAAAN1KKCyClhMR654y83XQYAAC0eF6sCAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwJqAgkpOTo169eqlt27Zq27at+vXrp7feest/3ufzafLkyYqNjVWbNm1066236uuvv7a9aAAA0DwEFEQ6dOigefPmaevWrdq6dauuv/563Xbbbdq1a5ckaerUqVq9erVyc3O1ceNGVVdX65ZbblFtba0jxQMAgOAWGsjgESNG1NmfM2eOcnJytHnzZnXo0EGLFi3SK6+8oqFDh0qSli5dKo/Ho3Xr1mnYsGH2VQ0AAJqFRl8jUltbq9zcXNXU1Khfv34qLCzU8ePHlZaW5h+TmJio5ORkbdq06azz+Hw+eb3eOhsAAGgZAg4iO3bs0MUXXyy32617771Xq1evVvfu3VVWVqbw8HC1a9euzvi4uDiVlZWddb7s7GxFRUX5N4/HE3gXAAAgKAUcRLp27aqioiJt3rxZ9913nzIyMrR79+6zjrcsSy6X66zns7KyVFlZ6d9KSkoCLQkAAASpgK4RkaTw8HBdfvnlkqTU1FRt2bJFTz/9tEaPHq1jx47p0KFDdV4VKS8vV//+/c86n9vtltvtbkTpAAAg2J33OiKWZcnn8yklJUVhYWHKy8vznystLdXOnTt/NIgAAICWK6BXRGbOnKn09HR5PB5VVVUpNzdX+fn5Wrt2raKiojR+/HhNmzZNMTExio6O1vTp09WzZ0//p2gAAAB+KKAg8t1332ncuHEqLS1VVFSUevXqpbVr1+rGG2+UJD355JMKDQ3VqFGjdOTIEd1www166aWXFBIS4kjxAAAguAUURBYtWvSj5yMiIrRgwQItWLDgvIoCAAAtA981AwAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMCSiIZGdn66qrrlJkZKTat2+vkSNHau/evXXGDB48WC6Xq842ZswYW4sGAADNQ0BBpKCgQBMnTtTmzZuVl5enEydOKC0tTTU1NXXGTZgwQaWlpf7tueees7VoAADQPIQGMnjt2rV19hcvXqz27dursLBQ1157rf9469atFR8f36A5fT6ffD6ff9/r9QZSEgAACGLndY1IZWWlJCk6OrrO8WXLlik2NlY9evTQ9OnTVVVVddY5srOzFRUV5d88Hs/5lAQAAIJIQK+I/JBlWcrMzNTAgQOVnJzsP37XXXcpKSlJ8fHx2rlzp7KysrR9+3bl5eXVO09WVpYyMzP9+16vlzACAEAL0eggMmnSJH366afauHFjneMTJkzw/zk5OVldunRRamqqPvnkE/Xt2/eMedxut9xud2PLAAAAQaxRb81MnjxZa9as0fr169WhQ4cfHdu3b1+FhYVp3759jSoQAAA0XwG9ImJZliZPnqzVq1crPz9fSUlJ57zNrl27dPz4cSUkJDS6SAAA0DwFFEQmTpyo5cuX6/XXX1dkZKTKysokSVFRUWrVqpW++OILLVu2TMOHD1dsbKx2796tadOmqU+fPhowYIAjDQAAgOAV0FszOTk5qqys1ODBg5WQkODfVqxYIUkKDw/Xu+++q2HDhqlr166aMmWK0tLStG7dOoWEhDjSAAAACF4BvzXzYzwejwoKCs6rIAAA0HLwXTMAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwJiAgkh2drauuuoqRUZGqn379ho5cqT27t1bZ4zP59PkyZMVGxurNm3a6NZbb9XXX39ta9EAAKB5CCiIFBQUaOLEidq8ebPy8vJ04sQJpaWlqaamxj9m6tSpWr16tXJzc7Vx40ZVV1frlltuUW1tre3FAwCA4BYayOC1a9fW2V+8eLHat2+vwsJCXXvttaqsrNSiRYv0yiuvaOjQoZKkpUuXyuPxaN26dRo2bNgZc/p8Pvl8Pv++1+ttTB8AACAIndc1IpWVlZKk6OhoSVJhYaGOHz+utLQ0/5jExEQlJydr06ZN9c6RnZ2tqKgo/+bxeM6nJAAAEEQaHUQsy1JmZqYGDhyo5ORkSVJZWZnCw8PVrl27OmPj4uJUVlZW7zxZWVmqrKz0byUlJY0tCQAABJmA3pr5oUmTJunTTz/Vxo0bzznWsiy5XK56z7ndbrnd7saWAQAAglijXhGZPHmy1qxZo/Xr16tDhw7+4/Hx8Tp27JgOHTpUZ3x5ebni4uLOr1IAANDsBBRELMvSpEmTtGrVKr333ntKSkqqcz4lJUVhYWHKy8vzHystLdXOnTvVv39/eyoGAADNRkBvzUycOFHLly/X66+/rsjISP91H1FRUWrVqpWioqI0fvx4TZs2TTExMYqOjtb06dPVs2dP/6doAAAATgkoiOTk5EiSBg8eXOf44sWLdffdd0uSnnzySYWGhmrUqFE6cuSIbrjhBr300ksKCQmxpWAAANB8BBRELMs655iIiAgtWLBACxYsaHRRAACgZeC7ZgAAgDEEEQAAYAxBBAAAGEMQAQAAxhBEAACAMQQRAABgDEEEAAAYQxABAADGEEQAAIAxBBEAAGAMQQQAABhDEAEAAMYQRAAAgDEEEQAAYAxBBAAAGEMQAQAAxhBEAACAMQQRAABgDEEEAAAYQxABAADGEEQAAIAxBBEAAGAMQQQAABhDEAEAAMYQRAAAgDEEEQAAYAxBBAAAGEMQAQAAxhBEAACAMQQRAABgDEEEAAAYQxABAADGEEQAAIAxBBEAAGAMQQQAABhDEAEAAMYQRAAAgDEEEQAAYAxBBAAAGEMQAQAAxgQcRDZs2KARI0YoMTFRLpdLr732Wp3zd999t1wuV53tmmuusa1gAADQfAQcRGpqatS7d28tXLjwrGNuuukmlZaW+rd///vf51UkAABonkIDvUF6errS09N/dIzb7VZ8fHyjiwIAAC2DI9eI5Ofnq3379rriiis0YcIElZeXn3Wsz+eT1+utswEAgJbB9iCSnp6uZcuW6b333tPjjz+uLVu26Prrr5fP56t3fHZ2tqKiovybx+OxuyQAAHCBCvitmXMZPXq0/8/JyclKTU1Vp06d9Oabb+r2228/Y3xWVpYyMzP9+16vlzACAEALYXsQOV1CQoI6deqkffv21Xve7XbL7XY7XQYAALgAOb6OyMGDB1VSUqKEhASn7woAAASZgF8Rqa6u1ueff+7f379/v4qKihQdHa3o6GjNnj1bv/zlL5WQkKCvvvpKM2fOVGxsrH7xi1/YWjgAAAh+AQeRrVu3asiQIf79U9d3ZGRkKCcnRzt27NDLL7+s77//XgkJCRoyZIhWrFihyMhI+6oGAADNQsBBZPDgwbIs66zn33777fMqCAAAtBx81wwAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwJtR0Ac3Nnj17bJ8zNjZWHTt2tH1eAABMI4jYpLb6kORyaezYsbbPHdGqtfb+Zw9hBADQ7BBEbHLSVy1ZlmJumaawGI9t8x4/WKKDbzyuiooKgggAoNkhiNgsLMYjd/zlpssAACAocLEqAAAwhiACAACMIYgAAABjCCIAAMCYgIPIhg0bNGLECCUmJsrlcum1116rc96yLM2ePVuJiYlq1aqVBg8erF27dtlWMAAAaD4CDiI1NTXq3bu3Fi5cWO/5+fPn64knntDChQu1ZcsWxcfH68Ybb1RVVdV5FwsAAJqXgD++m56ervT09HrPWZalp556SrNmzdLtt98uSVqyZIni4uK0fPly3XPPPedXLQAAaFZsXUdk//79KisrU1pamv+Y2+3Wddddp02bNtUbRHw+n3w+n3/f6/XaWVKz4cTS8RLLxwMAzLI1iJSVlUmS4uLi6hyPi4vTgQMH6r1Ndna2Hn74YTvLaFacXDpeYvl4AIBZjqys6nK56uxblnXGsVOysrKUmZnp3/d6vfJ47FsiPdg5tXS8xPLxAADzbA0i8fHxkv73ykhCQoL/eHl5+RmvkpzidrvldrvtLKNZYul4AEBzZOs6IklJSYqPj1deXp7/2LFjx1RQUKD+/fvbeVcAAKAZCPgVkerqan3++ef+/f3796uoqEjR0dHq2LGjpk6dqrlz56pLly7q0qWL5s6dq9atW+vOO++0tXAAABD8Ag4iW7du1ZAhQ/z7p67vyMjI0EsvvaQHHnhAR44c0e9//3sdOnRIV199td555x1FRkbaVzUAAGgWAg4igwcPlmVZZz3vcrk0e/ZszZ49+3zqAgAALQDfNQMAAIwhiAAAAGMcWUcEwcWJVVtZsRUA0BAEkRbMyVVbWbEVANAQBJEWzKlVW1mxFQDQUAQRsGorAMAYLlYFAADGEEQAAIAxBBEAAGAMQQQAABhDEAEAAMYQRAAAgDF8fBdoAsXFxaqoqLB9XlawBRDsCCKAw4qLi9W125U6euSw7XOzgi2AYEcQARxWUVGho0cOs4ItANSDIAI0EVawBYAzcbEqAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjCGIAAAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADAm1HQBaL727NnjyLyxsbHq2LGjI3MDAJoWQQS2q60+JLlcGjt2rCPzR7Rqrb3/2UMYAYBmgCAC2530VUuWpZhbpiksxmPr3McPlujgG4+roqKCIAIAzQBBBI4Ji/HIHX+56TIAABcwLlYFAADGEEQAAIAxBBEAAGAMQQQAABhjexCZPXu2XC5XnS0+Pt7uuwEAAM2AI5+a6dGjh9atW+ffDwkJceJuAABAkHMkiISGhvIqCAAAOCdHgsi+ffuUmJgot9utq6++WnPnzlXnzp3rHevz+eTz+fz7Xq/XiZKAcyouLlZFRYXt8zq11D0ANAe2B5Grr75aL7/8sq644gp99913+utf/6r+/ftr165diomJOWN8dna2Hn74YbvLAAJSXFysrt2u1NEjh02XAgAtiu1BJD093f/nnj17ql+/frrsssu0ZMkSZWZmnjE+KyurznGv1yuPx95lwYFzqaio0NEjhx1Zlv7Il1tV+f5SW+cEgObC8SXe27Rpo549e2rfvn31nne73XK73U6XATSIE8vSHz9YYut8ANCcOL6OiM/n0549e5SQkOD0XQEAgCBjexCZPn26CgoKtH//fn300Uf61a9+Ja/Xq4yMDLvvCgAABDnb35r5+uuvdccdd6iiokKXXHKJrrnmGm3evFmdOnWy+64AAECQsz2I5Obm2j0lAABopviuGQAAYAxBBAAAGEMQAQAAxhBEAACAMQQRAABgDEEEAAAYQxABAADGEEQAAIAxBBEAAGAMQQQAABhDEAEAAMYQRAAAgDG2f+kd0BT27NlzQc/XlJyqPTY2Vh07dnRkbgA4hSCCoFJbfUhyuTR27FjTpRjn9M8iolVr7f3PHsIIAEcRRBBUTvqqJctSzC3TFBbjsW3eI19uVeX7S22bryk49bOQpOMHS3TwjcdVUVFBEAHgKIIIglJYjEfu+Mttm+/4wRLb5mpqdv8sAKApcbEqAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGNYRwTAWTmxfLzP55Pb7bZ9Xoll6ZtKcXGxKioqHJnbqccwGGtuKQgiAM7g6PLxrosk66T984pl6ZtCcXGxuna7UkePHHZkficew2CsuSUhiAA4g9NL6bMsffCqqKjQ0SOHg+oxDMaaWxKCCICzcmopfZalD37B+BgGY80tARerAgAAYwgiAADAGIIIAAAwhiACAACMIYgAAABjCCIAAMAYPr4LAA3g1MqcTq7K6UTNTqy2i/q1lNVgCSIAcA5Orszp1KqcTq8mCme1pNVgCSIAcA5Orczp5KqcTtV8anVcOKslrQZLEAGABgrGlTmdWh0XTSMY/84FiotVAQCAMQQRAABgDEEEAAAYQxABAADGOBZEnnnmGSUlJSkiIkIpKSl6//33nborAAAQpBwJIitWrNDUqVM1a9Ysbdu2TYMGDVJ6erqKi4uduDsAABCkHPn47hNPPKHx48frd7/7nSTpqaee0ttvv62cnBxlZ2fXGevz+eTz+fz7lZWVkiSv12t7XdXV1f+7z7LPdfLYUVvnPvWRNrvndmpeJ+em5qaZm5pPm/u/X0uSCgsL/f/W7bJ3715JDvw8grHmIHwMnfpZSMFdc3V1ta2/a0/NZVlWYDe0bObz+ayQkBBr1apVdY5PmTLFuvbaa88Y/9BDD1mS2NjY2NjY2JrBVlJSElBusP0VkYqKCtXW1iouLq7O8bi4OJWVlZ0xPisrS5mZmf79kydP6r///a9iYmLkcrlsrc3r9crj8aikpERt27a1de4LUUvrV2p5PdNv89bS+pVaXs/NqV/LslRVVaXExMSAbufYyqqnhwjLsuoNFm63W263u86xn/zkJ06VJUlq27Zt0D/ggWhp/Uotr2f6bd5aWr9Sy+u5ufQbFRUV8G1sv1g1NjZWISEhZ7z6UV5efsarJAAAoGWzPYiEh4crJSVFeXl5dY7n5eWpf//+dt8dAAAIYiGzZ8+ebfekbdu21Z///GddeumlioiI0Ny5c7V+/XotXrzY8bddziUkJESDBw9WaGjL+L6/ltav1PJ6pt/mraX1K7W8nltav6dzWVagn7NpmGeeeUbz589XaWmpkpOT9eSTT+raa6914q4AAECQciyIAAAAnAvfNQMAAIwhiAAAAGMIIgAAwBiCCAAAMCZog8iGDRs0YsQIJSYmyuVy6bXXXjvnbQoKCpSSkqKIiAh17txZzz777BljnnnmGSUlJSkiIkIpKSl6//33nSg/YE70m52drauuukqRkZFq3769Ro4c6f+ipQuBU4/xKdnZ2XK5XJo6daqdZTeaU/1+8803Gjt2rGJiYtS6dWv97Gc/U2FhoRMtBMSJfk+cOKE//elPSkpKUqtWrdS5c2c98sgjOnnypFNtBCTQnktLS3XnnXeqa9euuuiii876d3XlypXq3r273G63unfvrtWrVztRfsCc6PeFF17QoEGD1K5dO7Vr105Dhw7Vxx9/7FQLAXHq8T0lNzdXLpdLI0eOtLNs44I2iNTU1Kh3795auHBhg8bv379fw4cP16BBg7Rt2zbNnDlTU6ZM0cqVK/1jVqxYoalTp2rWrFnatm2bBg0apPT0dBUXFzvVRoM50W9BQYEmTpyozZs3Ky8vTydOnFBaWppqamqcaiMgTvR8ypYtW/T888+rV69edpfdaE70e+jQIQ0YMEBhYWF66623tHv3bj3++OPG1/ORnOn30Ucf1bPPPquFCxdqz549mj9/vh577DEtWLDAqTYCEmjPPp9Pl1xyiWbNmqXevXvXO+bDDz/U6NGjNW7cOG3fvl3jxo3TqFGj9NFHH9lZeqM40W9+fr7uuOMOrV+/Xh9++KE6duyotLQ0ffPNN3aW3ihO9HvKgQMHNH36dA0aNMiOUi8sjfyS3QuKJGv16tU/OuaBBx6wunXrVufYPffcY11zzTX+/Z///OfWvffeW2dMt27drBkzZthXrA3s6vd05eXlliSroKDAljrtZGfPVVVVVpcuXay8vDzruuuus+6//37b6z1fdvX74IMPWgMHDnSkRjvZ1e/NN99s/fa3v60z5vbbb7fGjh1rX7E2aUjPP3S2v6ujRo2ybrrppjrHhg0bZo0ZM+a8a7STXf2e7sSJE1ZkZKS1ZMmS8ynPdnb2e+LECWvAgAHWiy++aGVkZFi33XabXWVeEIL2FZFAffjhh0pLS6tzbNiwYdq6dauOHz+uY8eOqbCw8IwxaWlp2rRpU1OWaotz9VufyspKSVJ0dLTj9TmhoT1PnDhRN998s4YOHdrUJdqqIf2uWbNGqamp+vWvf6327durT58+euGFF0yUe94a0u/AgQP17rvv6rPPPpMkbd++XRs3btTw4cObvN6mcrafSzA+bzXG4cOHdfz48aB93mqIRx55RJdcconGjx9vuhRHtJj1ZMvKys740r24uDidOHFCFRUVsixLtbW19Y45/Qv8gsG5+k1ISKhzzrIsZWZmauDAgUpOTm7KUm3TkJ5zc3P1ySefaMuWLYaqtE9D+v3yyy+Vk5OjzMxMzZw5Ux9//LGmTJkit9ut3/zmN4Yqb5yG9Pvggw+qsrJS3bp1U0hIiGprazVnzhzdcccdhqp23tl+LsH4vNUYM2bM0KWXXhr0/7E4mw8++ECLFi1SUVGR6VIc02KCiCS5XK46+9b/Lyrrcrnq/Pn0MacfCxY/1u/pJk2apE8//VQbN25sktqc8mM9l5SU6P7779c777yjiIgIE+XZ7lyP8cmTJ5Wamqq5c+dKkvr06aNdu3YpJycn6IKIdO5+V6xYoaVLl2r58uXq0aOHioqKNHXqVCUmJiojI6PJ620qzel5KxDz58/XP/7xD+Xn5zebf9M/VFVVpbFjx+qFF15QbGys6XIc02KCSHx8/Bn/QygvL1doaKhiYmJkWZZCQkLqHXP6/zaCwbn6/aHJkydrzZo12rBhgzp06NCUZdrqXD2/+eabKi8vV0pKiv98bW2tNmzYoIULF8rn8ykkJKSpy260hjzGCQkJ6t69e50xV155Zb0X8F7oGtLvH//4R82YMUNjxoyRJPXs2VMHDhxQdnZ2sw0iZ/u5BOPzViD+9re/ae7cuVq3bt0FddG5nb744gt99dWo9LJ0AAADKUlEQVRXGjFihP/YqU+AhYaGau/evbrssstMlWebFnONSL9+/ZSXl1fn2DvvvKPU1FSFhYUpPDxcKSkpZ4zJy8tT//79m7JUW5yrX+l//2uaNGmSVq1apffee09JSUkmSrXNuXq+4YYbtGPHDhUVFfm31NRU3XXXXSoqKgqqECI17DEeMGDAGR/J/uyzz9SpU6cmq9MuDen38OHDuuiiuk9rISEhF8zHd51wtp9LMD5vNdRjjz2mv/zlL1q7dq1SU1NNl+OYbt26nfGcdeutt2rIkCEqKiqSx+MxXaI9zFwje/6qqqqsbdu2Wdu2bbMkWU888YS1bds268CBA5ZlWdaMGTOscePG+cd/+eWXVuvWra0//OEP1u7du61FixZZYWFh1quvvuofk5uba4WFhVmLFi2ydu/ebU2dOtVq06aN9dVXXzV5f6dzot/77rvPioqKsvLz863S0lL/dvjw4Sbvrz5O9Hy6C+lTM070+/HHH1uhoaHWnDlzrH379lnLli2zWrdubS1durTJ+zudE/1mZGRYl156qfXGG29Y+/fvt1atWmXFxsZaDzzwQJP3V59Ae7Ysyz8+JSXFuvPOO61t27ZZu3bt8p//4IMPrJCQEGvevHnWnj17rHnz5lmhoaHW5s2bm7S3+jjR76OPPmqFh4dbr776ap3nraqqqibtrT5O9Hu65vipmaANIuvXr7cknbFlZGRYlvW/B+u6666rc5v8/HyrT58+Vnh4uPXTn/7UysnJOWPev//971anTp2s8PBwq2/fvhfMR1md6Le++SRZixcvbpqmzsGpx/iHLqQg4lS///rXv6zk5GTL7XZb3bp1s55//vkm6ObcnOjX6/Va999/v9WxY0crIiLC6ty5szVr1izL5/M1UVc/rjE91ze+U6dOdcb885//tLp27WqFhYVZ3bp1s1auXNk0DZ2DE/126tSp3jEPPfRQk/V1Nk49vj/UHIOIy7L+/2ovAACAJtZirhEBAAAXHoIIAAAwhiACAACMIYgAAABjCCIAAMAYgggAADCGIAIAAIwhiAAAAGMIIgAAwBiCCAAAMIYgAgAAjPk/Gfrk4GUwUc8AAAAASUVORK5CYII=",
"text/plain": [
"Figure(PyObject )"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"┌ Warning: `getindex(o::PyObject, s::Symbol)` is deprecated in favor of dot overloading (`getproperty`) so elements should now be accessed as e.g. `o.s` instead of `o[:s]`.\n",
"│ caller = top-level scope at In[16]:6\n",
"└ @ Core In[16]:6\n"
]
}
],
"source": [
"# \n",
"# A COMPLETER \n",
"# "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 10\n",
"\n",
"Adapt the code such that each edge has now a cost function of the form\n",
"$\\ell_e(x_e) = t_e(1+0.15*(x_e/c_e)^4)$.\n",
"\n",
" Plot an histogram of the price of anarchy over $100$ realization of $t_e$ and $c_e$. Comment."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# \n",
"# A COMPLETER \n",
"# "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## More complex graphs\n",
"\n",
"### Question 11\n",
"\n",
"Construct a graph with $6$ nodes. Define linear costs on each edge. \n",
"Compute the price of anarchy."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# \n",
"# A COMPLETER \n",
"# "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 12 \n",
"\n",
"On your new graph test if it is possible to add an edge with null cost\n",
"that increase the cost of the user equilibrium."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# \n",
"# A COMPLETER \n",
"# "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Question 13\n",
"\n",
"By introducing edge-flux variables $x^1$ and $x^2$ rewrite\n",
"the user equilibrium and social optimum problems in the case\n",
"where $K=2$. \n",
"\n",
"Construct an example on your previous graph and compute the price\n",
"of anarchy."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# \n",
"# A COMPLETER \n",
"# "
]
}
],
"metadata": {
"@webio": {
"lastCommId": null,
"lastKernelId": null
},
"kernelspec": {
"display_name": "Julia 1.4.0",
"language": "julia",
"name": "julia-1.4"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.4.0"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": true,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": true,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 2
}