diff --git a/Pollutant_Problem.ipynb b/Pollutant_Problem.ipynb new file mode 100644 index 0000000..a0e7e03 --- /dev/null +++ b/Pollutant_Problem.ipynb @@ -0,0 +1,37 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "hello world" + ] + } + ], + "metadata": { + "anaconda-cloud": {}, + "kernelspec": { + "display_name": "Python [default]", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +} diff --git a/Population.ipynb b/Population.ipynb new file mode 100644 index 0000000..884795f --- /dev/null +++ b/Population.ipynb @@ -0,0 +1,339 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Partial Differential Equation Uses in Aging Population Predition" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The ageing population in any countries is a challenge. Many countries in the world are facing this unavoidable social problem. Ageing population occurs when the median age of a country grows caused by the higher life expectancy and decline in the birth rate or/and fertility rate. \n", + "According to data from World Population Prospects: the 2015 Revision (UN, 2015), the total number of those aged 60 years or over persons who are defined as the older person has grown recently in most countries, and this development is projected to be higher in the coming years.\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Modeling population" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To implement appropriate policy to decrease ageing population, the current situation of this part of people should be consider. However, as we can imagine, any people at any age may go on to next age safely or die. Also some women at certain age can give birth to a baby. So, setting up a appropriate mathematical model about ageing population seem to be complicated. \n", + "\n", + "Here we use the McKendrick PDE to simulate the population model:\n", + "\n", + "$$ \\frac{\\partial p(t,x)}{\\partial t} + \\frac{\\partial p(t,x)}{\\partial x} = -d(t,x)p(t,a) $$\n", + "\n", + "$$ t = t_0 : p = p_0(x) $$\n", + "\n", + "$$ x = 0 : p(t,0) = b(t) $$\n", + "\n", + "In this partial differential equation, $p(t,x)$ represents the distribution of at time $t$ and age $x$, $d(t,x)$ is the fertility rate of those people whose age is $x$, and $b(t)$ is the birth rate in a ragion at time $t$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Problem set up" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let’s use this model to simulate the population development in the U.S. from 2006 to 2016. We assume that the highest age in this model is 99. According to the data from United Nations Population Division, we assume that the death rate is 0.008, and the birth rate is 0.013\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 44, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import numpy \n", + "from matplotlib import pyplot\n", + "from matplotlib import rcParams\n", + "rcParams['font.family'] = 'serif'\n", + "rcParams['font.size'] = 16" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here, we use time-forward space-backward scheme to solve the problem. Then we can re-write the equation as: \n", + "\n", + "\\begin{equation}\\frac{p_i^{n+1}-p_i^n}{\\Delta t} + \\frac{p_i^n - p_{i-1}^n}{\\Delta x} = -dp_i^n, \\end{equation}\n", + "\n", + "We solve for this unknown as follows:\n", + "\n", + "$$ p_i^{n+1} = (1 - d\\Delta t - \\frac{\\Delta t}{\\Delta x}) p_i^n + \\frac{\\Delta t}{\\Delta x} p_i^{n-1} $$\n", + "\n", + "In this model, we have boundary condition which is: \n", + "\n", + "$$ x = 0 : p(t,0) = b(t) $$\n", + "\n", + "\n", + "\n", + "Because we choose 2006 as the first year of the model, so the initial condition will be the population distribution in the U.S.by age in 2006(numbers in thousands).\n", + "\n", + "Here, we divide the population into 13 group according to the chart. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "![population_2006](population2006b.png)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 73, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "nx = 101 \n", + "dx = 100/(nx-1)\n", + "\n", + "nt = 20 \n", + "dt = 0.5 \n", + "\n", + "d = 0.008\n", + "b = 0.013\n", + "\n", + "x = numpy.linspace(0,100,nx)" + ] + }, + { + "cell_type": "code", + "execution_count": 74, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "#initial condition\n", + "p = numpy.ones(nx)\n", + "\n", + "p[:4] = 20363/5\n", + "p[5:9] = 19629/5\n", + "p[10:14] = 20651/5\n", + "p[15:19] = 20916/5\n", + "p[20:24] = 20393/5\n", + "p[25:29] = 20138/5\n", + "p[30:34] = 19343/5\n", + "p[35:44] = 43121/10\n", + "p[45:54] = 42797/10\n", + "p[55:64] = 30981/10\n", + "p[65:74] = 16554/10\n", + "p[75:84] = 12962/10\n", + "p[85:] = 3989/15" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Boundary condition\n", + "\n", + "In this problem, we should consider the boundary condition. That is, in every time step, the amount of people whose age is zero--called new baby--is equal to the number of total in last time step times the birth rate.\n", + "\n", + "Because we assume the birth rate is a constant, so the boundary is a dirichlet boundary condition:\n", + "\n", + "$$ x = 0 : p(t,0) = bp_{total}(t-1) $$" + ] + }, + { + "cell_type": "code", + "execution_count": 76, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "def population(p, b, d):\n", + " \"\"\"Computes populatoin density\n", + "\n", + " Parameters\n", + " ----------\n", + " p : popualtion density of the first year\n", + " b : birth rate\n", + " d : death rate \n", + " Returns\n", + " -------\n", + " p: \n", + " popualtion density of the last year\n", + " \"\"\"\n", + " for n in range(1, nt): \n", + " pn = p.copy() \n", + " p[2:] = (1 - d - dt/dx)*pn[2:] + (dt/dx)*pn[1:-1] \n", + " \n", + " amount = numpy.sum(p)\n", + " p[0] = amount*b\n", + " \n", + " return p " + ] + }, + { + "cell_type": "code", + "execution_count": 77, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "p_2016 = population(p, b, d)" + ] + }, + { + "cell_type": "code", + "execution_count": 78, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAELCAYAAADkyZC4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4FFX3wPHvSaH3XqWDIBCaiiAQkCIgRYpiQURFfEWx\nCyoqyKvIi10sPxQQCxYQG0oVQ1Ok9yaIFBGQ3iHl/P7Yia7rJtmE3WzJ+TxPnmXunpk5M8CezNx7\nZ0VVMcYYY/wtKtgJGGOMiUxWYIwxxgSEFRhjjDEBYQXGGGNMQFiBMcYYExBWYIwxxgSEFRhjjDEB\n4VOBEZFKInJCRFa6/axyXgu5xeUXkbEisllE1ovITBGp42V7MSIyUkQ2ichaEVkkIs3T2Pf9IrJB\nRFaLyHIR6Zb1wzXGGJNdYjIRu0xV22QQMxXID8Sp6jkReQZIEJE4Vf3DLW4sEA80U9XDInI7MFtE\nrlDVtalBIjIUeBC4TFV/E5G2wHci0kVVZ2Uid2OMMdlMfJnJLyKVgInpFRgRaQfMAlqr6nynLRbY\nB0xW1XudtprAJuA2VZ3ktv56YIeqdnGWCwO/A2NUdYRb3HSgkqrWy+zBGmOMyT7+7IPpCZwHFqc2\nqGqis9zTLa6H85rgsf48oL2I5HOWOwJ504ir4xQqY4wxISozBaaMiHwgIj87fSwfiUhdt/frAXtV\nNcljvR1AaREp4RaXAuzyEhcD1HGLS233jAOon4ncjTHGZDNfC0wykAi8qKqXA02c5Z9FpLETUwI4\n4WXd485rcbe40/rve3Pe4vCyzeOAuMUZY4wJQT4VGFXdo6pxqrraWT4J3AWcAp4LYH7GGGPCVJb7\nYFT1LLAOaOo0HQQKeglNHcZ8yC0un4iID3F42aZnnDHGmBDk0zBlZ67LGafT3l0yEO38eS3QWERi\nPPphqgD7VfWgW1wfoCL/7IepAiQBG93iACp7iVO39z1ztS+4McaYLFBVz1/8L4ivVzCv8s+RYKlD\nkOsBK5ymaUAs0Mwjpjmu+TGpvnBe4z320RqYpaqnneWZwBkvcW2Ajaq6Na1kVdV+VHn66aeDnkOo\n/Ni5sHNh5yL9n0DIzC2yh0WkDICIRAMv4OqIHwGgqnNwzYMZKSJ5nXWG4boqGZW6EXUVhnHAYyJS\n3Nlef6Aq8IRb3DFgJDBIRKo4cW2BdsBDmT5SY4wx2crXmfwvAAOBGU7fSXFckyWvUtUFbnG9gNHA\nahFJAvYA8frPWfwA9wBPA4tF5DyukWLtVHWde5CqjhaRM8B0EUnEdUuul6rOztRRGmOMyXY+FRhV\n3QAM9iHuNHCvD3HJwFPOT0axrwGv+ZCm8RAfHx/sFEKGnYu/2bn4m52LwPLpUTHhREQ00o7JGGMC\nTUTQIHXyG2OMMZliBcYYY0xAWIExxhgTEFZgjDHGBIQVGGOMMQFhBcYYY0xAWIExxhgTEFZgjDHG\nBIQVGGOMMQFhBcYYY0xAWIExxhgTEFZgjDHGBIQVGGOMMQFhBcYYY0xAWIExxhgTEFZgjDHGBIQV\nGGOMMQFhBcYYY0xAWIExxhgTEFZgjDHGBIQVGGOMMQFhBcYYY0xAWIExxhgTEFZgjDHGBIQVGGOM\nMQFhBcYYY0xAWIExxhgTEFZgjDHGBIQVGGOMMQFhBcYYY0xAZKnAiMhCEUkRkYv8nZAxxpjIkOkC\nIyI9geaAenkvv4iMFZHNIrJeRGaKSB0vcTEiMlJENonIWhFZJCLN09jf/SKyQURWi8hyEemW2ZyN\nMcZkv0wVGBGJBUYB36YRMhWoD8Spal1gKZAgImU94sYCvYHmqlofmAjMFpH6HvsbCjwOdFbVBsBQ\nYIqIdMhM3sYYY7KfqP7rQiTtYJEHgMbAL8BTQBVV3eW81w6YBbRW1flOWyywD5isqvc6bTWBTcBt\nqjrJbdvrgR2q2sVZLgz8DoxR1RFucdOBSqpaL40cNTPHZIwxBkQEVRV/btPnKxgRKQY8DDyWRkhP\n4DywOLVBVROd5Z5ucT2c1wSP9ecB7UUkn7PcEcibRlwdp1Cl6Zdd+1mwcmt6IcYYYwIoM7fIngLe\nV9XdabxfD9irqkke7TuA0iJSwi0uBdjlJS4GqOMWl9ruGQeuW3Fp+mjGUlrf9RJvfPZDemHGGGMC\nxKcCIyI1gF7Ac+mElQBOeGk/7rwWd4s77eU+lrc4vGzzOCBucV7tP3ycB25sy9jPErj/xU9JTk5J\nL9wYY4yf+XoF8zwwSlW9FZCQtP/wcZrWrcJPE4ewfvteuj/8JidOnQ12WsYYk2PEZBQgIi2AusB1\n7s1eQg8CnqPFAAo5r4fc4vLJv3vjvcUBFASOpBP3L8OHD2f53MUUPr6REjE9mfHaYAaNnsyVd4zh\nm5cHcVGZYmmtaowxOUJCQgIJCQkB3UeGo8hEZARwI3/fqhKgDFAK12iw87iGEncD+gMF3PthRORr\noImqlnOWh+C61fbXCDSn/TVgIFBUVU+LyPXAZFyj0ha4xT0IjAFqq+q/evFT61aNa59k+suDqFW5\nDACqyiuTv+eFD+cwbcxdXF63iu9nyRhjIlxQRpGp6tOqWkNVGzk/DYG3nbc7Om0zgWlALNDMLeFY\nXJMyp7pt8gvnNd5jV62BWap62lmeCZzxEtcG2OituLg7cPg4pYsX+mtZRHjgpra8/diNXHP/WD6e\nuTTd4zbGGHNhsvosstQq99f6qjoH1zyYkSKS12keBiThmpyZGrcVGAc8JiLFAUSkP1AVeMIt7hgw\nEhgkIlWcuLZAO+Ch9JI7ey6Rs+eTKFwg77/e69Iyju/feoDH3/ySYW9+SUqKdf4bY0wgZHYmf0cR\nWQXc6TR9KyIr3UJ6AWuB1SKyAWgKxKvqHx6bugeYAiwWkbXAHUA7VV3nHqSqo4FngekishoYDfRS\n1dnp5bn/8HFKFSuIiPervfo1KvDze0OZv/IXej76f5w8bZ3/xhjjb5mayR8ORESXrt/Bf57/iOUf\nPJFu7Lnzidz9/Mcs37STL1/4D1XKl0g33hhjIlVQZ/KHk/2Hj1O6WKEM43LniuXdJ/tyW9dmXHHb\naH5YviUbsjPGmJwhMgvMIdctMl+ICPfdcBUfPnMbfR5/hzc++4FIu6ozxphgiMgCc+DICZ+uYNy1\nvbw2P04YwltTF3Dnsx9y7nxigLIzxpicISILjK+3yDxVq1CSnyYO4dCxk7S+6yX+OHgsANkZY0zO\nEJkF5tAJn2+ReSqYPw9TRw+kY7O6XNZvFMs2/Obf5IwxJoeIyAJz4EjWrmBSRUVF8eQdnRn7SB86\n3fc6n8xa5sfsjDEmZ8jwWWThaP+hE/+YxZ9V3eIbULlccbo99BabfvuDpwdcQ1RURNZkY4zxu4j8\ntNx/+DilimbtFpmnuJoV+fm9ocz5eRM3PzmBxKRkv2zXGGMiXUQWmKMnTlOiSAG/ba908ULMe+tB\njp86S+8h/2cjzIwxxgcRWWCKFspPdLR/Dy1P7limjbmLmOhouj30JqfPnvfr9o0xJtJEZIEpncUR\nZBnJFRvDJ8/dQYkiBeh83+t2JWOMMemIzALjhw7+tMTERDNpeH+KFMzHQy9PzXgFkyMlJSXz+4Ej\n9lQIk6NFZIHxVwd/WqKjo3hv+K3MWrKRyfa9Msax7+Axrrl/LNW6PUH+FoOp3Xs4g8d8YkXG5FgR\nOUz5QubA+KpwgbxMHX0nbe9+hbgaFbikWrmA79OErpSUFPoNf4/aVcrw0gO9qVS2GGfOJdLqzhd5\ndvx3DLujc7bl8vHMpcxduvmv5UvrVOKuXq2ybf/GpIrIK5hA3iJzF1ezIi/c34uej77NiVP2nTI5\n2asfz+P4qTO8cF8valYqTe5csRQpmI+Zrw9m4jc/Mm7agow34gezl2zkwZencHndyjSrX5Ur6lfl\n+Umz+GbBmmzZvzHuIvIKJtC3yNz1u+YK5i3bzH/Hf8vowT2zbb8mdKzespvnJs7g5/eGEhMT/Y/3\nypYozKyx99FywAuULlaIbvENApbHtt0H6PvUBKaOHkiLhjX+aq9duQw9h/wfKy6+iPKligZs/8Z4\nsisYP3huUHfe/Woxvx84kq37NcF3+ux5bnjiXV5+sDdVK5T0GlO9Yik+ee4O7nvxs4BN1D1+8gxd\nH3yTZwZ2/UdxAWjeoDr39I7npmETSE62rwg32ScyC0yAhimnpXypotzRrTnPvPNttu7XBN+Y92cR\nV6MCN3dqmm5cy0Y1qVq+BB/N+NnvOagq/Ya/R6tGNRjYs6XXmMf6d0QEnps4w+/7NyYtEVlgSmVD\nJ7+nIf2uZtoPq9i6c3+279sER2JSMuO+WMSw2zv5FD/s9k6MmjjT71cR81dsZcOve3n14evTjImO\njuKj/97Oa5/MY8fvB/26f2PSEpkFJhv7YFIVK5yfB29qy5Nvf5Xt+zbBMX3hWqqWL0Hd6uV9im/d\npBbFi+Rn6vcr/JrH85NmMuSWDuSKTb9LtVzJIvRp34QPvlvi1/0bk5aILDB5cscGZb+D+7Rh4apt\nrNi0Myj7N9nr7c8XcFcat6S8ERGeuK0Tz06YQUqKf65iVm7exfrte7m50+U+xfe75gre/3aJzc0x\n2SIiC0yw5M+bm2G3d2LEuOnBTsUE2PY9f7Jqy256tmmUqfU6Na9LdFQU0xeu80seoyfN5MEb25I7\nl2+/VDWuXYlcsTH8uGa7X/ZvTHqswPjZrV2asXjtdnbvOxzsVEwAjZu2gH7XNM301bLrKqYjz074\n7oKvIn7ZtZ95y7Yw4NoWmdp/v85NmfTtTxe0b2N8YQXGz/LlycUN7S9l/FeLg52KCZBz5xOZ+M1P\n3JmJD3Z3Pdo05NCxU/y8fscF5THmg9nc3bsVBfPnydR6N3e6nKnfr+SMPRHcBJgVmAC4s0cLxn+9\nmCT7crKI9Pm8VcTVqECNi0pnaf2oqCju6tmSN6fMz3IOe/88ytTvV3Lv9W0yvW75UkVpUrsSX9vs\nfhNgVmACoH6NClQoVYSZP20IdiomAMZNW5ipzn1v+ndpxjcL1/LnkRNZWn/sZz9w09WXZfmL9fpd\ncwWTpttoMhNYVmAC5M5rWzBu2sJgp2H8bM/+I6zb/jvXtKh3QdspXqQA3ePjmJCFW6mnz57nnS8W\nMbhP5q9eUl3buiFL1v/KvoPHsrwNYzJiBSZArmvXhEVrtrFnvz0+JpJ8Nmc53Vs18HnUVnoG9Y7n\nrc8XZHri5YffLeGKelWzfIsOXH2F3Vs14NM5y7O8DWMyYgUmQPLnzc0NHS5lwtfW2R9JPp2znD7t\nm/hlW03qVKZU0YLM+HG9z+uoKq98PI8Hbmp7wfu/pkU9u41rAsoKTADdeW0L3v1ykT1gMEL8uudP\nduw9ROsmtfy2zUHXxWeqs3/2ko3ExkQT37jmBe/7qstqs3jNdhtNZgLGCkwAxdWsSKliBZm3bHPG\nwSbkfTpnOb2uavSvR/JfiOvaNmbZxt/YtvuAT/GvfPw999/QBhG54H0XLpCX+tXLs3D1tgveljHe\nZFhgRKSqiLwgIstFZJmIbBGRBSLSySMuv4iMFZHNIrJeRGaKSB0v24sRkZEisklE1orIIhFpnsa+\n7xeRDSKy2tl/t6wfanDc0tn1aA4T/j6Z7b/bY6ny5snFvde35ok3v8wwdtOOP1i1ZTc3dLjMb/vv\ncMUlzLLbZCZAfLmC6QhcB/RW1UtVtRawCPhaRNxnmk0F6gNxqloXWAokiEhZj+2NBXoDzVW1PjAR\nmC0i9d2DRGQo8DjQWVUbAEOBKSLSIdNHGUR92jfhm4VrOXnavvEynG3a8QcHj57kygbV/b7th/u2\nZ8m6HSxYuTXduNGTZjGwRwu/PmuvQ9M6zPppo9+2Z4w7XwrM78BwVXWfdjzaWbcbgIi0AzoAT6rq\nOSdmJBCNq0jgxNUEBgCjVPUwgKqOB3YAz7rFFQaGAW+o6m9O3FxgNvBCpo8yiEoVK0TLhjWYNm9V\nsFMxF+DT2cu5rm1joqL8f1c5X55c/G9wD+574bM0++u+WbCG+Su38uBN7fy678a1K7Hv0DEb7WgC\nIsP/Lar6papO8Ggu7Lym3jjuCZwHFrutl+gsu3+PcA/nNcFje/OA9iKSz1nuCORNI66OU6jCxi2d\nm/K+PSI9bKkqn8xeRp8OlwZsH9e1a0KBfLm9jjrcd/AYdz77IR88cxuFC+T1636jo6Noe1ltZi+x\nqxjjf5n+dUxEyuO6zbUceNNprgfsVdUkj/AdQGkRKeEWlwLs8hIXA9Rxi0tt94wD1624sHFNi/qs\n2rLbHoAZplZu3sX5pGQuu6RywPYhIrz68PU8+fbXHDt55q92VaX/M5O4o/uVAbk9B9DhijrMWmL9\nMMb/fC4wTmf/L7iKQxRwraqedN4uAXh75sVx57W4W9xp/fdjZL3F4WWbxwFxiwsLeXLH0vuqxnw0\nc2mwUzFZ8H/TFnB71+Z+GbmVnkYXX0SXFvXp9uCbvDJ5LotWb+PFD+dw+NgpnhpwTcD2275pHeYu\n3WzD6Y3f+VxgVPVXVa2B6/bYL8BaEWkWsMwiTN9Ol9sXPYWh4yfPMGXuSm7r5nWgo9+9/GBv+l3T\nlF92H+DBl6bw+qc/8NF/byfWj0OjPZUvVZRyJQqzbONvAduHyZnS/45VL5yrlgec0VxvAg2Ag4Dn\naDGAQs7rIef1IJBPRMTjKsZbHEBB4Eg6cWGjWVw1zp1PZPnGnVwawFstxr8mz1xKmya1KFuicMbB\nflAgXx76d21O/67ZU9BSpQ5Xblqvarbu10S2DAuMiORRVW9jbNcBPUUkFlgLNBaRGI9+mCrAflVN\nLRhrgT5ARf7ZD1MFSAI2usUBVPYSp27vezV8+PC//hwfH098fHx64dlCRBhwbQvGfvYDk0b0D3Y6\nxgeqyv9NW8jowT0yDg5zHZrW4elx3/D0nV2CnYrJJgkJCSQkJAR2J6qa7g/wA3C5l/alwGHnz+2A\nZKCl2/uxuK40XnNrq+nE3eKxrXXA127LhYGTwFMecdOBdRnkq6Hq8LGTWrT1/bp73+Fgp2J8sHT9\nDq3a9XFNTk4OdioBd+bseS3YcrAePnYy2KmYIHE+OzOsCZn58bUPZoSIFEtdEJHBQGPgVecTfQ4w\nCxgpIqnjKIfhuioZ5VbMtgLjgMdEpLizrf5AVeAJt7hjuObRDBKRKk5cW6eQPeRjziGnaKH83NK5\nKa99Mi/YqRgf/N+0BQy4tkVA5r6Emjy5Y7kyrjpzl9pjjYz/+PI/53Fcky0TRGSliGzCNRP/JlUd\n4RbXC9etq9UisgFoCsSr6h8e27sHmAIsFpG1wB1AO1Vd5x6kqqNxTb6cLiKrcU3u7KWqszN9lCHk\ngRvbMv7rxRx3G4pqQs+xk2f4fN4q+nfJOeNYrm52CTN/tOHKxn9EI2xU07/HD4SeG594l0YXX8TD\nfdsHO5UMnT57noWrfiFhxVYuKlOM9k3rUK1CyWCnFXAvfTiHJet38NnzdwY7lWyzded+Wt/1Enu+\nez7gQ7JN6BERVNWvf/GRf+0fgh7u255XP5nH+UTPeamhITEpmU9mLaP9oFco3f5h/jv+O2Jjovl5\n/Q5a3DGG6t2H8cQbX3L0xOlgpxoQ23YfYNR7MxkxMGd1eNe4qBS5c8WwYfveYKdiIkSmhymbC9fo\n4ouoeVEpPp29nL6dmwZ0Xys37+LJt75i+aadxERHEx0lFMqfl+Zx1WjdpBYtG9UgV2wMR0+c5sjx\n08xespG3Pp9PjYqluLt3PFNHD6SQ2+NJVJV1237ntU/mUbPHUwzt14G7e8dn+ABGVeV8YpJfvgky\nkJKSkrnl6Yk8eUdnalfxNvI+cokIV19xCTN/2kDd6uWDnY6JAHaLLEjmLdvMrcPfY/H4R6lYpljG\nK2TSlt/2Meytr1i8ZjvDbu/Eta0bkpKSQlJyCoeOnWLByq0krNjKotXbUKBIgbwULZSPBjUrcu/1\nrYmrWTHDfWz8dS+Pjf2SZRt/o+1ltWnRsDrN6lfj9NnzrNv2O+u2/c6WnfvZte8wO/cd4tz5JC6t\nU5nu8XF0j29wQV/5GyjPjv+OH1ZsYfbY+3JE576nrxJW8/pnPzD3zQeCnYrJZoG4RWYFJohe/HAO\n479azMJ3HqZ4kQJ+2WZKSgovfTSX0ZNm8dDN7bj3+tbkz5vbL9tOy5bf9jF/5VYWrtrGj2u3Uyh/\nXupVL0+96uW4uHIZKpUtTqWyxcmTK4aEFVv54ofVfDV/NRdXLsPTA64h3o/fEHkhVm7exdX3vsaK\nDx4PSNEPBydOnaVcx0fZN2tMwP/dmNBiBcYH4VRgAB599XMWrv6FuW8+cMH/oXftO0y/pyeSlJzC\nB8/0p3K5EhmvFCSJScl8NONnnp0wg3IlCzOodzz1qpenesVSxMZEc/jYKX5a9ys/rtnOlp372b3/\nMLv3HyEmOoqbO13ObV2bU71iKb/kkpKSwufzVvHwK1MZdc+13Hi1/77QKxy1HvgiD/dtT+cr62Uc\nbCKGFRgfhFuBUVX6j5jE3j+Pct8NV1GrUmkqly2OiLD/8HF27TvM0ROnKV64ACWLFqBk0YL/KETJ\nySksWr2Nz+et5JPZy3nwprY80rc90dHhcXsnKSmZj2ct47O5y9n8m6uQlCxakGMnz3DZJZVpHleN\nS6qWo2LpolQsXYyjJ04z8Zsf+eC7n7mkaln+N7hnlh+9o6rM/XkTj73xBarw/D3X0q7pv76ENccZ\n/d5M9hw4wuuP3hDsVEw2sgLjg3ArMOD6bf6Zd6azbONvbNm5n32HjpOSohQrlI+KpYtRpGBeDh07\nxcGjJzlw5ATRUVGUKV6I0sUKsW3PAcqVKELPNg25vn2TkOzXyIyz5xL5/c+jVCpTjJh0HvB4PjGJ\nyTOXMnTsF9zQ/lJG/qcrBfLlyXD7W3fuZ/aSjSSs2MqCVb9QokgBht95Db2uapQj+1y8Wb1lN72H\njuOXL0YGOxWTjazA+CAcC4ynM2fPExUlXkdcqSonTp1l36Hj7Dt0jPKliuaIeSlpOXj0JI+8OpV5\ny7bwSN/2dG0Vx0Vu/SfJySls3PEHXyasZsrcFRw8epJOzevSqlFNWjWu+Y9Y46KqVOg0lLlv3p/j\nRtLlZFZgfBAJBcZk3oKVWxn/1WK+XbSOi8oUI65mBTb++gcbfv2DMsUL0al5XXq3bUzzuGp2peKD\nx8Z+wfnEJF58oHewUzHZxAqMD6zA5GxJScn8uHY7m3bs45Jq5ahfvfw/5vEY3/y6508uv/V5dn/7\nfIZznExksALjAyswxvhH+0GvcGuXZjl+VF1OYY+KMcZkmzuvbcG4aQuDnYYJY1ZgjDFedYtvwOad\n+9jy275Mr6uq7N53mMSk5ABkZsKF3SIzxqQpK539Kzfv4sGXprB22x5Onz1P5bLFqVO1LKPv7RH2\nw+gjmfXB+MAKjDH+s33PnzT1sbN/38FjDB37BbN+2sCIgV24rWtzkpJT2Lb7AF8vWMPEb37ip4lD\nKOGnxyIZ/7I+GGNMtqpWoSQNa1Vk0vSf0o1bsu5XmtzyHCWLFmDL589wZ4+WxMREkyd3LHWrl+fx\n2zrR66pGdHvwTc6eS8ym7E2w2RWMMSZdq7fspuPg1/jf4J5ev15i4teLGfL6F4x/si9dWsaluZ2U\nlBRufGI8AJOfvd3mI4UYu0XmAyswxvjfxl/30uHe1xjarwODrmuNqrJi007e/nwBC1b9wlcv3u3T\nrP+z5xK56j8v07VlfYbcenU2ZG58ZQXGB1ZgjAmMHb8fpN2gV4irUYEVm3eRKzaaXlc14pG+7Sla\nKL/P29n82z5a3fkCO78ZZZM4Q4gVGB9YgTEmcP44eIxPZi2j7eW1qVutHCJZ+zzqOPg1el/VmNu6\nNfdzhiarrMD4wAqMMaFv1k8bePS1aayePCzLRcr4l40iM8ZEhPZN65CYlEzCiq3BTsUEkBUYY0y2\nExHu69OGVyZ/H+xUTABZgTHGBEXfzk1ZvGYb2/f8GexUTIBYgTHGBEW+PLm4o/uVvP7pvGCnYgLE\nCowxJmgG9Y5n0vQlnDtvs/sjkRUYY0zQVCxTjIsrl7bO/ghlBcYYE1TdWzXgq/lrgp2GCQArMMaY\noOoW34CvF6zB5q9FHiswxpigurhyGfLnyc2KTTuDnYrxMyswxpig6x4fx5cJq4OdhvEzKzDGmKDr\nZv0wESnDAiMicSIyTkQ2isgaEVkvIq+KSAmPuPwiMlZENjsxM0WkjpftxYjISBHZJCJrRWSRiHh9\n4p2I3C8iG0RktYgsF5FuWT9UY0yourxuFQ4cOcGvNukyovhyBfMpUBRopKpxQDugPbBIRHK7xU0F\n6gNxqloXWAokiIjnl0SMBXoDzVW1PjARmC0i9d2DRGQo8DjQWVUbAEOBKSLSIbMHaYwJbdHRUXRt\nGWdXMRHGlwKTAjyqqmcBVPUPYAxQA+gEICLtgA7Ak6p6zllvJBCNq0jgxNUEBgCjVPWws73xwA7g\nWbe4wsAw4A1V/c2JmwvMBl7I4rEaY0JYt1ZxfDnf+mEiiS8Fpr6q7vBo2wsIrisbgJ7AeWBxaoCq\nJjrLPd3W6+G8Jnhsbx7QXkTyOcsdgbxpxNVxCpUxJoJcdenFrN6ym4NHTwY7FeMnGRYYVU3y0lwL\n15XNfGe5HrDXS+wOoLRbf009Z71dXuJigDpucantnnHguhVnjIkgefPkos2lFzNj8fpgp2L8JNOj\nyEQkCrgNeFdVtzvNJYATXsKPO6/F3eJOe/lGMG9xeNnmcVxXTsUxxkScjs0uYeZPG4KdhvGTrAxT\nfgrX7bAH/JyLMSaHu7pZXWYv2UhyckqwUzF+kKkCIyL9gV7A1ap6xu2tg0BBL6sUcl4PucXlk39/\nR6q3OLxs0zPOGBNBLipTjFLFCtqs/ggR42ugiPTFddXSWlU9P+DXAo1FJMajH6YKsF9VD7rF9QEq\n8s9+mCpAErDRLQ6gspc4dXvfq+HDh//15/j4eOLj49MLN8aEkI7N6jLjx/VcVrdKsFOJaAkJCSQk\nJAR0H+LLA+ZE5GbgUeAqVf3TaesMlFPVd5xhyjNxFZ8FzvuxwD7gI1Ud7LTVBDYB/VX1fbftrwN2\nqGpXZ7n9Ao9kAAASdUlEQVQw8DvwP1V9xi1uOlBJVVMHAXjL1UsXjzEmXMz9eRNPvv0VP00cGuxU\nchQRQVU97y5dEF9m8t8EjAPewzWU+CanrQtQFkBV5wCzgJEiktdZdRiuq5JRqdtS1a3Oth4TkeLO\n9vsDVYEn3OKO4ZpHM0hEqjhxbXFN8nzoAo7XGBPiWjSszsZf/+CQDVcOe77cInsNyI1rcqWnEW5/\n7gWMBlaLSBKwB4h3Jma6uwd4GlgsIudxjRRrp6rr3INUdbSInAGmi0gikAz0UtXZPuRsjAlTuXPF\n0qpRTeb8vIk+HS4NdjrmAvh0iyyc2C0yY8Lfm1MS+Hn9DiaN6B/sVHKMoNwiM8aY7Hb1Fa75MCkp\nNlw5nFmBMcaEnKoVSlKkQD5Wb90T7FTMBbACY4wJSR2bX2KPjQlzVmCMMSGpS4v6fJGwKthpmAtg\nBcYYE5JaNarJrn1H7EvIwpgVGGNMSIqJiaZH64ZMmbsi2KmYLLICY4wJWde1a8xnVmDClhUYY0zI\natmwBnsOHGG73SYLS1ZgjDEhKyYmmp5tGtltsjBlBcYYE9J6t23MZ3OWBzsNkwVWYIwxIa1lwxrs\nPXiMX3btD3YqJpOswBhjQlp0dBQ929hosnBkBcYYE/Kua9uET+eswB5kG16swBhjQt6VDapzPjGJ\nmT9uCHYqJhOswBhjQl50dBTPDerOY298YU9YDiNWYIwxYaF7fAPy5Irlk9k2oixcWIExxoQFEeH5\ne67lybe+4nxiUrDTMT6wAmOMCRvxTWpRs1Jpxk1bGOxUjA+swBhjwsqoQdfy7ITvOHHqbLBTMRmw\nAmOMCSsNalWkZ5tGtLrzBXtGWYizAmOMCTuvP9qH/l2acUX/0Xz+/cpgp2PSIJE2cUlENNKOyRjj\n3bINv3HdY+OoVak0FUoVpXjh/BTMl4cjJ06z//BxDhw+QYF8ualctjiVy5Wgad0qXFa3SrDTDkki\ngqqKX7cZaR/GVmCMyVmOnjjNvGWbOXTsFIeOneL4qTMUK5SfUkULUqpYQU6cPsdvew+y84/DTF+0\njvrVy/PcoO7UrV4+2KmHFCswPrACY4xJy7nzibw5ZT6j3ptJlxb1ef3RPuTLkyvYaYUEKzA+sAJj\njMnIsZNnuPv5yfxx8BjfvDyI/HlzBzuloLMC4wMrMMYYXyQnp3DbM5PYte8w01+5J8cXmUAUGBtF\nZozJkaKjo5jwVD8qlS3ONfeP5dSZc8FOKeLYFYwxJkdLTk7h1uHvcS4xiU9HDUDEr7/Ehw27gjHG\nGD+Ljo7inWF9+fX3P3ll8vfBTieiWIExxuR4eXLHMnX0QJ6fNJOFq34JdjoRwwqMMcYAlcuVYNLw\nW+nz+Lv8cfBYsNOJCD4XGBEpKyIzRcS+7ccYE5GublaXAd2v5Lqh4+wrAfzApwIjIj2AH4GqQJo9\n6CKSX0TGishmEVnvFKQ6XuJiRGSkiGwSkbUiskhEmqexzftFZIOIrBaR5SLSzcdjM8aYTHtqQGeK\nFcrPvf/7BBswdGF8vYJ5BGgLLM4gbipQH4hT1brAUiBBRMp6xI0FegPNVbU+MBGYLSL13YNEZCjw\nONBZVRsAQ4EpItLBx7yNMSZToqKi+OCZ/ixas423ps4PdjphzadhyiISpaopIjIRuEVVo73EtANm\nAa1Vdb7TFgvsAyar6r1OW01gE3Cbqk5yW389sENVuzjLhYHfgTGqOsItbjpQSVXrpZGrDVM2xlyw\n7Xv+pNlto/n0uQHEN6kV7HQCLhDDlGN8CVJVX/pdegLncbvKUdVEEVnsvHev09zDeU3wWH8eMFBE\n8qnqaaAjkDeNuDEiUlNVt/qSvzHGZFa1CiWZ/N/b6fPEu0z73100i6vm87opKSms3rqH75duYu7S\nzazfvpdihfJRsmhByhQvxH96taJFwxoBzD40+FRgfFQP2Kuqnj1jO4DOIlJCVQ86cSnALi9xMUAd\nYLkTl9ruGQeuW3FWYIwxAXPVZbUZ/+QtXPvIWzzevyOD+7RJdyLmoaMnmfD1Yt6cOp9cMTG0u7w2\n/+nVioa1KnLs5BkOHD7Btt0HuGnYeJrVr8aY+3pSsUyxbDyi7OXPAlMCOOGl/bjzWhw46MSd9nIf\nyz0udXt42eZxQNzijDEmYDpfWY8lE4fSa8j/8ePa7Yy5rxcVSxf9q9CcOHWWuUs38WXCar5esJau\nLevz6XMD0vzembaX1+aWa65g9KSZNLjpv7x4fy9u7dIsOw8p2/izwBhjTESqUr4Ei8c/ysOvTOXy\nfqM4dfY8tSuXIW/uXKzYvJMr6lWl85X1GHNfT0oVK5Th9vLlycWIgV25scNltBr4IhVKFaXt5bWz\n4Uiylz8LzEHAc7QYQOrZPuQWl0/+3RvvLQ6gIHAknbh/GT58+F9/jo+PJz4+PoPUjTEmfXlyxzJ2\nyA2MHXIDR46fYtOOfRw7eYYWDatTIF+eLG2zVuUyfDZqAL2HjmPBuIepVbmMn7NOW0JCAgkJCQHd\nR6YedpnBKLK3gP5AAfd+GBH5GmiiquWc5SHAc0AVVd3lFvcaMBAoqqqnReR6YDKuUWkL3OIeBMYA\ntb118tsoMmNMuBn/5SJGvz+LJROHUqxw/qDkEOoPu5wGxAJ/3Ux0hik3xzU/JtUXzmu8x/qtgVnO\nCDKAmcAZL3FtgI02gswYEylu734lXVvG0efxdyJqcmdmC0ya1U1V5+CaBzNSRPI6zcOAJGCUW9xW\nYBzwmIgUBxCR/rieEvCEW9wxYCQwSESqOHFtgXbAQ5nM2xhjQtroe3uw//AJvvhhVbBT8RtfJ1r+\nD9cHe0WgKLDGeesyj9th+YDRQHtchWUPcL+qbvLYXjTwNHAdrrkzJ4BHVPVHL/sejOvWWSKQDAxX\n1W/SydVukRljwtKcJRu5e/THbPjsaXLFZu8YLPvKZB9YgTHGhLOOg1+jY7O6DO7TJlv3G+p9MMYY\nYy7QmME9eXbCdxw9cTrj4BBnBcYYY0JI3erl6dYqjucmzAh2KhfMCowxxoSYEQO7Mv7rxez8I83p\nfmHBCowxxoSYsiUKM6D7lbz44Zxgp3JBrMAYY0wIGtynDR/O+JmDR08GO5UsswJjjDEhqFzJIvRo\n3ZA3Pvsh2KlkmRUYY4wJUY/0bc8bU+Zz+uz5YKeSJVZgjDEmRNWqXIbmcdWY+HVG31YfmqzAGGNM\nCHv0lva88OEckpKSg51KplmBMcaYEHZF/WpULF2UKXNXBDuVTLMCY4wxIW5Iv6sZ/f6ssHvSshUY\nY4wJcZ2a10UVZixeH+xUMsUKjDHGhDgR4bFbr+bZCTPC6irGCowxxoSB3m0bc+DICRau+iXYqfjM\nCowxxoSB6OgohvTrwHMTw+chmFZgjDEmTPTtdDnrt+9lxaadwU7FJ1ZgjDEmTOTOFcvDN7dj1MSZ\nwU7FJ1ZgjDEmjAy4tgVL1v/KotXbgp1KhqzAGGNMGMmfNzcv3t+Lu5+fHPKz+63AGGNMmLmuXRNK\nFyvE65+G9pOWrcAYY0yYERHGPtqHZyd8x94/jwY7nTRZgTHGmDBUq3IZ7uzRgodenhrsVNJkBcYY\nY8LUsNs789O6X/lk1rJgp+JVTLATMMYYkzX58uTi65fupv09r5I3dyzd4hsEO6V/sAJjjDFhrH6N\nCnz7yj10HPw6uXPFcHWzusFO6S92i8wYY8Jc49qV+PKF/9D3qYl8s2BNsNP5i4TTkzl9ISIaacdk\njDG+WLjqF24f+T7lSxZhxMAutGxU0+d1RQRVFX/mYwXGGGMiSFJSMh/NXMoz70ynfKkixDeuRYOa\nFYirWZFyJQqTJ3csIv+uI1ZgfGAFxhhjIDEpmW8XrWPFpp2s3rqbNVv38OfRk5xPTKJA3tzkzZML\nAaKiooiKEvZ8N9oKTEaswBhjTNoSk5I5deYcZ84loqqoKikpykVli1uByYgVGGOMybxA3CIL+VFk\nIlJSRD4Ukc0isklEpohI+WDnZYwxJn0hXWBEJBaYC8QCtYE6wCngBxHJF8zcjDHGpC+kCwxwK1AX\neFQdwBCgKvCfYCZmjDEmfSHdByMiM4CLVbWKR/ta4KSqNvOyjvXBGGNMJuXEPpj6wA4v7TuAetmc\nS9hJSEgIdgohw87F3+xc/M3ORWCFeoEpAZzw0n4cyCciubM5n7Bi/3n+Zufib3Yu/mbnIrBCvcAY\nY4wJU6FeYA4CBb20FwJOq+q5bM7HGGOMj8Khk7+Wqlb1aE+3kz+78jPGmEji707+UP8+mGnA2yJy\nkaruAhCR0rjmxAzxtoK/T5AxxpisCfUrmFhgGbAJuBlQYDzQDGioqqeDmJ4xxph0hHQfjKomAu2A\nZGAjsAEoALSx4mKMMaEtpK9gjMkKEVkINAcqp95aNSYnEZGywESgvaoG7UIipK9gfJUTH4gpInEi\nMk5ENorIGhFZLyKvikgJj7j8IjLWOTfrRWSmiNQJVt6BJiI9cRWXf/3mlJPOhYj0FJH5IrJMRLY7\nrze5vZ8jzoWINBGR70Rkg/P/5GcR6eURE1HnQkR6AD/ieqRWmlcQvh63iMSIyEjns3WtiCwSkeY+\nJZP6fQDh+oPrQZhrgE8BcX7eA7YC+YKdXwCPezMwBcjjLJfF1Ve1GcjtFjcDWJDaBjwDHADKBvsY\nAvRvYSvwDa7bqhd5vJ8jzgXwALA89biAaOBDYEJOOhdAJeCo83mQerdmIJACdI7UcwH8BFTDdQWT\nnE6cT8cNvO18rhRzlm/H9dDh+hnmEuyT4YeTOcD5MKnk1lYaSAIeCnZ+ATzujUAVj7bbnHNxrbPc\nzvnP1MotJhY4BLwe7GMIwDl5wPkgfdqzwOSUcwFUBs4CjTzay6S25aBz8R/n30F9j/ajwEeRei6A\nKOc1zQLj63EDNZ1z2M9j/fXANxnlEgm3yHoAu1R1Z2qDqu7H9QHcM2hZBV59VfV8TtteXFdwRZ3l\nnsB5YHFqgLoGTiwmws6NiBQDHgYeSyMkp5yLvsARVV3p3qiq+9zacsq5SHJeYz3ao/i7eyDizoWq\npvgQ5utx93BeEzzWnwe0z+hrUyKhwOTIB2KqapKX5lq4fiuZ7yzXA/Z6id0BlPbsrwlzTwHvq+ru\nNN7PKefiCuA3EekhIgucPrrFItLfLSannItPcN3aGeb0N4iIPAHkwnXbB3LOufDk63HXw/WZ4jlY\nZgeueZTp9lVFQoGxB2ICIhKF6xbZu6q63WlO79wAFM+O3AJNRGoAvYDn0gnLEecCqIjrO5QeAnqq\nah3gZWCciKRe3eWIc6GqJ4C2QF5cj53aj+s7ptqpauovYTniXHjh63GXwPVYLs/BAj6dn0goMMbl\nKVyXvA8EO5EgeB4Y5Xyg5HR5gHzAw6r6J4CqTgW+Ah4XkbzBTC47iUhNYCmu37aLqGopYBjwhYh0\nCGpyOUQkFJgc/0BM5/ZHL+BqVT3j9lZ65wZcHXphTURa4PqN/W33Zi+hEX8uHKlFdo1H+ypchac2\nOedc/BcoDNyf+jmgqp/iGjk1ybnqzynnwpOvx30Q150gz/9TPp2fSCgwa3GNnPFUBViXvalkPxHp\ni+uqpbWqev5lrwXKiYjnM+eqAPtV9WB25BhgbXH9O14mIitFZBWuoagA3zltV5MzzgW4+hzg3/+3\nk93ac8q5qAvs8fJL5lagJK7jzSnnwpOvx70W17+Zil7iknANpkpTJBSYaUAlEbkotcHtgZhTg5ZV\nNhCRm4FHgKtSb4eISGcRGeCETMM1gqaZ2zqxuCYiRsS5UdWnVbWGqjZyfhry99VMR6dtJjngXDi+\ncV7re7TXA87getxSTjkXB4CyzpWKu8q4JiAeIeecC0++HvcXzmu8x/qtgVma0SO7gj1m2w9jvmOB\n1cDHuCaUReEa/72FyJ5oeRNwGnjQ+XPqz9vAU25x3+EaVZbXWR6Bq7MzLCeR+XhuhuMxNyqnnAvn\n3/8SXMNI8zttLXDNjRmaw85FT+ffwX/d2loD54API/1c4Jpgmt5ES5+OG3gL1yTu4s5yf1wTLetl\nlENEPItMREriGilzKa4hdetx3Xf9PaiJBZCIHAKKpPH2CFV9xonLB4wG2uO6pN2D69xsypZEs5GI\ndMQ1kqy087MJOK+qjZz3c8S5EJEiuI6zA66rlnPAa6o6wS0mp5yLdsBQXBNNk3F9PryPazJhohMT\nUedCRP6HayJlRVxz4lL74y5Tt2HJvh63iETjmrx8Ha6BRCeAR1T1xwxziYQCY4wxJvREQh+MMcaY\nEGQFxhhjTEBYgTHGGBMQVmCMMcYEhBUYY4wxAWEFxhhjTEBYgTHGGBMQVmCMMcYEhBUYY4wxAWEF\nxhhjTED8P7xsEaCESGr0AAAAAElFTkSuQmCC\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "pyplot.plot(x, p_2016, color='#003366', ls='-', lw=1)\n", + "pyplot.ylim(0,5000);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Comparsion \n", + "\n", + "Let's compare this result to the actual population density in 2016\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 90, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "z = numpy.linspace(0,100,101)\n", + "p_2016a = numpy.ones(101)\n", + "\n", + "p_2016a[:4] = 20201/5\n", + "p_2016a[5:9] = 16348/5\n", + "p_2016a[10:14] = 17677/5\n", + "p_2016a[15:19] = 17040/5\n", + "p_2016a[20:24] = 16585/5\n", + "p_2016a[25:29] = 18101/5\n", + "p_2016a[30:34] = 18962/5\n", + "p_2016a[35:44] = 40071/10\n", + "p_2016a[45:54] = 43007/10\n", + "p_2016a[55:64] = 36482/10\n", + "p_2016a[65:74] = 21713/10\n", + "p_2016a[75:84] = 13061/10\n", + "p_2016a[85:] = 7493/15" + ] + }, + { + "cell_type": "code", + "execution_count": 91, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAELCAYAAADkyZC4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXecFOX9+N+fLVc4qoCKWMCaIII1KsR4FiyxRUGjSYyx\nJv7UxC7GBhI1JqZoxG8ssSUmGo0lagQ05ixoLFFERYIaFBFFkHLAtS3P74+ZuZvbnb3bKbszt/u8\nX697HTf77DOfGXafz3zqI0opNBqNRqMJmljYAmg0Go2mMtEKRqPRaDQlQSsYjUaj0ZQErWA0Go1G\nUxK0gtFoNBpNSdAKRqPRaDQlQSsYjUaj0ZSEohSMiGwlIutE5A3bz5vm74G2cQ0icrOILBSRd0Rk\nloiMcZgvISIzROQ9EZkvIi+KyMQC5z5XRN4VkXki8rqIHOX9cjUajUZTLhIuxr6mlNq/lzEPAQ3A\neKVUu4hcDTSJyHil1Ge2cTcDjcAEpdQqETkVmCMieyul5luDRGQqcD7wNaXURyJyIPAPETlCKTXb\nhewajUajKTNSTCW/iGwF3NWTghGRScBsYD+l1HPmsSTwOfBnpdQ55rHtgfeAU5RS99je/w6wWCl1\nhPn3IOBT4JdKqem2cU8AWymldnJ7sRqNRqMpH0HGYCYDHcBc64BSKmX+Pdk27hjzd1PO+58FDhKR\nfubfhwL1BcaNMRWVRqPRaCKKGwWzqYj8UUReMWMs94nIWNvrOwHLlFLpnPctBjYRkWG2cVlgicO4\nBDDGNs46njsOYJwL2TUajUZTZopVMBkgBfxKKbUnsLv59ysisps5ZhiwzuG9zebvobZxLSrfN+c0\nDoc5mwGxjdNoNBpNBClKwSilliqlxiul5pl/rwd+BGwAri2hfBqNRqPpo3iOwSil2oC3gb3MQyuB\nAQ5DrTTmL23j+omIFDEOhzlzx2k0Go0mghSVpmzWurSaQXs7GSBu/ns+sJuIJHLiMKOB5UqplbZx\nxwNb0D0OMxpIAwts4wBGOYxTttdzZdUb3Gg0Go0HlFK5D/6+KNaCuZHumWBWCvJOwH/MQw8DSWBC\nzpiJGPUxFo+YvxtzzrEfMFsp1WL+PQtodRi3P7BAKbWokLBKKf2jFFdddVXoMkTlR98LfS/0vej5\npxS4cZFdKCKbAohIHLgBIxA/HUAp9TRGHcwMEak333M5hlVynTWJMhTDbcClIjLUnO9kYGvgMtu4\ntcAM4CwRGW2OOxCYBFzg+ko1Go1GU1aKreS/Afgh8JQZOxmKUSx5gFLqedu4KcD1wDwRSQNLgUbV\nvYof4GzgKmCuiHRgZIpNUkq9bR+klLpeRFqBJ0QkheGSm6KUmuPqKjUajUZTdopSMEqpd4EfFzGu\nBTiniHEZ4Erzp7exNwE3FSGmJofGxsawRYgM+l50oe9FF/pelJaiWsX0JUREVdo1aTQaTakREVRI\nQX6NRqPRaFyhFYxGo9FoSoJWMBqNRqMpCVrBaDQajaYkaAWj0Wg0mpKgFYxGo9FoSoJWMBqNRqMp\nCVrBaDQajaYkaAWj0Wg0mpKgFYxGo9FoSoJWMBqNRqMpCVrBaDQajaYkaAWj0Wg0mpKgFYxGo9Fo\nSoJWMBqNRqMpCVrBaDQajaYkaAWj0Wg0mpKgFYxGo9FoSoJWMBqNRqMpCVrBaDQajaYkaAWj0Wg0\nmpKgFYxGo9FoSoJWMBqNRqMpCVrBaDQajaYkJMIWQKOJOus2tHHaz+6lI5Upy/m+1bgzJx2+d1nO\npdGUEq1gNJpe+HTFGl6c9yE3X3x8yc/14rwPeGruO1rBaCoCrWA0ml5IpzMMGdCPo/fbpeTnymSy\n3D/n9ZKfR6MpBzoGo9H0QiqdIZEoz1clkYiTSpfHFafRlBqtYDSaXkhnsiTiZVIw8RjpjFYwmspA\nKxiNphdS6QzJRLws50pqC0ZTQXhSMCLygohkRWTLoAXSaKJGOpMhES+PgjEsmGxZzqXRlBrXCkZE\nJgMTAeXwWoOI3CwiC0XkHRGZJSJjHMYlRGSGiLwnIvNF5EURmVjgfOeKyLsiMk9EXheRo9zKrNH4\nQVswGo03XCkYEUkC1wFPFhjyEDAOGK+UGgu8CjSJyIiccTcDxwITlVLjgLuAOSIyLud8U4GfAocp\npXYGpgIPisjBbuTWaPxQ3hhMXMdgNBWD22/N2RhKIy+PUkQmAQcDVyil2s3DM4A4hpKwxm0PnA5c\np5RaBaCU+gOwGLjGNm4QcDkwUyn1kTnuGWAOcINLuTUaz2gLRqPxRtEKRkQ2Ai4ELi0wZDLQAcy1\nDiilUubfk23jjjF/N+W8/1ngIBHpZ/59KFBfYNwYU1EV5P0ly3n+jUU9DdFoisKwYHQMRqNxixsL\n5krgXqXUJwVe3wlYppRK5xxfDGwiIsNs47LAEodxCWCMbZx1PHccGK44R5rXt3LnY3PZ70e/5ld/\nnEPz+tZuP+tb2gq9VaPJI11mCyad1gpGUxkUVckvItsBU4Cv9jBsGLDO4Xiz+XsosNIc16KUyk0S\nsI+z5sNhzmZAbOPy2PywqbS1p4jHYlx808NcOvNRamsSiAgALW0dvHL3VHb76lY9XI5GY5BKZ8oX\ng9EuMk0FUWyrmJ9jxEycFEjkaH7uRgCy2SzNG9o48co7+fnZx7DjNpsBcMCZv2bNupYwRdQUyT9e\nfJt/zH2nLOdKJuLMOPNI+ver63Y8ncmW14LRQX5NhdCrghGRfYCxwHH2ww5DVwK52WIAA83fX9rG\n9RMRybFinMYBDABW9zAuj2nTpnX+u7Gxkcd/c3a313Ugte/wl9mvAbDn2NElP9e1dz3FKUdNZKdt\nR3Y7XtZWMfGY/mxqykJTUxNNTU0lPUcxFsyBGLGa10wXkwCbmq/9Q0Q6MLLE5gO7iUgiJw4zGliu\nlLIUxnzgeGALusdhRgNpYIFtHMAoh3HK9noedgXjhJEKqv3cfYF0JssR3xjHdw75WsnPdcdjLzou\n7ulMmWMw+rOpKQONjY00NjZ2/j19+vTAz9HrY5lS6iql1HZKqV3Nn12A35svH2oemwU8DCSBCdZ7\nzbqZiRj1MRaPmL8bc061HzBbKWX5rmYBrQ7j9gcWKKU8p4hpC6bvUO4U4bTD58KIwZQvi0x/NjWV\ngle733KRdb5fKfU0MBuYISL15uHLMayS62zjFgG3AZeKyFAAETkZ2Bq4zDZuLUYdzVkiMtocdyAw\nCbjAo9yAsZC8OO8DLr/lUbJZ/bQYZaJQg1L+GIz+TGoqA7eV/IeKyJvAGeahJ0XkDduQKRiuq3ki\n8i6wF9ColPosZ6qzgQeBuSIyHzgNmKSUets+SCl1PUbx5RMiMg+4HpiilJrTk5zpdIa29lTB15OJ\nONtuMZzn3nifyRffqtOWI0wUFExZs8ji2rrWVA6uvjVKqaeUUrsopTZTSsWVUmOVUrvaXm9RSp2j\nlNpBKbWjUupgpdR7DvNklFJXKqW+opQap5SaqJR6qcA5bzLn2lkptZtS6vHe5Hxj4RK+ftovCr6e\niMfoV1fLM7ecy0YDG5h46i9Z/OnKguM14VHexd3ZPaXb9Ws03qjIdv3LVzWzyUYDC75uPanW1iS5\n44oTOeXICex9yvW88Ob7ZZRSUwxRCLBHwYrSaPoiFbll8her1hWlYABEhJ+ccAA7br0ZAxvqCr5H\nEw5RWNzL2a4/Ho+RzSqy2SyxWEU+/2mqiIpUMMtXNbPxRgMKvu60kBy4Z09NCjRhEQUFk0pnaKir\nLYsMItLZj6xGKxhNH6ciFcwXq9ex1aYFO8kYX2Dthig5T774NkecN7OosaccOYE7rvh+3vEoKJh0\nunwxGDDaxaQzWWqSZTulRlMSKlLBtLal2HRocS6y3vjk81VsselGQYlWVSz/spmTDt+LPzgoDjtP\nvvg2tz78guNrUcjgKmccCArX42g0fY2KtMFvvex7nNBD5XexCmbJ56vY9XvXcL/ZrkTjjlQ6Q00i\nQSwW6/GnJpko+P8RhT5gRquY8ikYXWypqRQqUsH0RjIRJ1VEKuiWm27EM7ecy9SbH+GqW/+uizJd\nUqx7qyeFHwkXWRmVnCWHLrbUVAJVqWDcFLON334LXrl7Kk+/8h7fu+JO/WTpgkpRMOV004EuttRU\nDlWpYNxu6rTJ0IE8+3/n07yhjfN+9dcSSlZZVIqCCceC0QpG0/epyCB/b3gpZqurTfLwL3/EitV9\nYkucSBCUggm70WT5LRhnOZ54YT7PvrawLDLU1iSZdsbh1OpUNo0PKlLBrFnXwuAB/Qq+7rVauiaZ\nYOTGQ/yI1mdQSvH0K+/12NPNYtjg/kwYv03e8SAUTBQaTYaSReYgx11/f4kBDXWMy9mvphRcc9dT\nnPatr7PN5sNLfi5N5VKRCmbyxb/nn/93fsHXdZZO7yxbsYYjz5/JQXuO6XFcOpPlxbc+6NxF1E4q\nnaFfXU2v50om4nSk0o6vRcFF5sWKSqczxOOxzm263dCTJXXMfrtw5L7jXc/pllsffkF/RzS+qUgF\n01ObGNBZOsXQ3pFm06GD+PtvzupxXGtbB4P3O8/xtY50mkGJesfX7EQ+BpMuzor6fOVaTvvZH3lv\n8Wcs/WINHzwyw1MNVcGeaGXuy6YVjMYvFRnkL0bBBPXlmXbr47z74bJA5ooSqXSGmmTvi1lPNSyV\nEuRPZXqPwWSzWU6adjfbbjGcp276Mc3P/dZzgW6hLLKOVLqoe/GXWa/y9L8X9DquJ3qyKjWaYqlM\nC6aHKn4IVsGMHjmMyRf/ntfu/SkDGupYtmINtzzYhFK9v3fPsaPL4u7wQrELe9xceDOZbOe/3c5R\n6P8jm82ilCIWc+9m8kLBdv3pbK+Fljf+5VnWtbRxw0+m+C7KTCScWxkVcz/n/HsB5//mQWbf/BNf\nMmgLRhMEFWnBbDykcKNLMHo9BfXlOenwvdlz7Gh+9ocnAXhx3gc82vQW/epqevz5bOVabnrg2UBk\nKAUdLiyHQk+7fhWM9X4vcQwv+GnXX5OM86erTwmk4r+n+9GTVfnBJ19w4pV38tefn8G47Tb3JUNN\nUisYjX8q0oLZfJOeM72CrjO49qxvMe6EGfz4+P1JpTOM335zLjv1mz2+57n/LOLKW/8emAxBY7V5\nKQZrMcqNtgSlYMpFz+36e34WO+u4/Xp8fUNrO6fN+CP3Xn1yr9dkdVPOpaf70by+lSPPv4Wrf3gk\n++yyXd7rK9es52d/eJJfnXtsnqXphLZgNEFQkRbMQXv1nPkU9Jdn5MZDOO2oiVx9+5NF+8mj7uN2\ns7j7VRDRVzD+U6Ub6mv5YnUzf571alFyOD0AdaQK34+Tp9/Dvrtuxw8nf8Px9SED+vHWoqVce9dT\nRckb9c+npm9QkQqmN0rxdHbJSYdQkzS+lMUFx6P9hFisogRzMapgBRNUoeVlp3yTa+98ikwvGYyF\ngvyGi8zZqrzu7KO58cJvF5wzHo9x389O5ab7ny1qe/Cekjc0mmKpSgVTijqYjQY18LuLTyj6addY\nzKKbKl1sFhkUXozcKJh0xgjo20lnyrwPS8F2/cEUe+63+w4MHdzA3559o8dxXmJB22+1SUHlY7HZ\n8MEcf9Du/PEf/+5VVu0i0wRBRcZgeqOUdTA9uTFyZSjkgrjlwSZmPthU1Pl+duZRHL3fLm5ELIpy\nushEhLgZd7CPD8OCKdiu36HQctZL77D/Hl/pdWG3EBEuO+Wb/HTmoxx74G4FkxcKPQC5sSoL8f3D\n9uaEy+7gitMO6zF5QisYTRBUrYIp1Zen2OB4Ty6ItxYtZcoBu3Lcgbv3OMdN9z/Lgv99VhIF48pF\nFi+cRVa0FWT+n4StYIrdcOzDpSv4/lV3s+SJ61yd45sTx/K7Bww31dYF2rD4yWbrjd3HbMXWI4ex\ndPnqHut0jAcgrWA0/tAKJmCCCPJ3pNOMGjGUHbfZrMc5Nhs+iI50aQKxPfn7cykUT/JrBUVFwTjF\nYG792/OcdPhe1NW6awYpIsz6Xc81Kl5iMG7OP2fmub2Oi3qMUNM30DGYgElljKf2tetbexzXkwXj\nzs1WOkusXC6yQnNERcHkuu7aO1Lc/cTLnHH0PiWTo5hCyzsfm8vS5atLJoNWMBq/VK0FU8oYTEcq\nxdhvT2fxY9cULLwrlHkFxT+pGkqqxfG1s6//C2/+95Ne50gkYvz5Z6fmdYmOgoIxgvzhb1WcSnff\nMvlvz77J+O02Z7stNymrHPbi12Ur1nDhjQ+xaN8ZJZFBKxhNEFStgillDGbEsMGMHD6YWS+/y+H7\njHMcV9ODDEaqc3EKppCbbdbL73LN/zuKzXvZXuCMa//EJ8tX5ymYYq0o8F/Jb80RBQvGsV1/jhy3\nPfwC53y758LKUshhvx83//VffPeQrzFscP+SyaDrYDR+0QrGRjqd4UfX3cf6lvZe5xg+pD+/u/iE\nvOPWInDG0ftw28MvFFQwPSm5IFxkHak0e4/bhi17abg4dFCD/wB9wTTl4tN7nRa0qLjIcmMwl5x0\nMPvvsUPJ5Egk8i0YpRQZM227pa2D2x95kZfuvLhkMtQkdB2Mxj9VqWAKBVHXrm/lgadf5/bLTuzx\n/UopvnP5H7jpouPzUj2tIP+3D9qdC298iKXLVzu2runpCdHIRCuiWLMnJVXsHMmEo5IKwkXmtlgz\nChZMMTGYQyeODeR8q5s38N3L7+Tx35zVrX2LU7q0vS/bH5/8N3vvtLVvF93sl99l7fpWjpuUn62o\nXWSaIKhKBVOwFUc6Q//6Wo4/eI9e5/j+VXc5FuBZ8ZOG+lqOP2gP7vz7XK48/fC898fjMZRy7kLc\nkfbvInNT5Og0h+tK/hIoKSP2UcZCyxJX8ucyZGADK9es56mX3ulm6To9ANn/P259+Hl+dd6xvs+/\nobWdux5/qaCCKRQj1GiKpSqzyHp64i4+Ndd5cbe7t86csm/BrZtFxPeTv7EIFEh1LjaOU8AV4qXZ\npdMcvoL8EbJgSpVscNZxjdzy4HPdjjk1u7Tfyzkzz6Vxt+19n/uAr32VuW99SGtbR95rOk1ZEwRa\nwdgIQsHY3Vs7bTuSHx+/fw9zeGvL3tv7wVB0xc7h2EfMxe6JJU1TLmMWmZEe3H1hz2aNv4vpQOyF\n4w7cjdcWfMQHn3zRXQ4HBWN9NocN7h/IFgaD+tczbtuRvDDvg7zXtItMEwS9fmtEZGsRuUFEXheR\n10TkvyLyvIh8M2dcg4jcLCILReQdEZklInltjUUkISIzROQ9EZkvIi+KyMQC5z5XRN4VkXnm+Y/y\nfqldFG7F4Saw7ewWcu9acraCineR5cuglCpajmIssd6opCyylEPsw3KPrW7eEPg56+tqOOfb+3HZ\nLY92HnP6fJbqXhy01xhmv/xu3nGtYDRBUMxj2aHAccCxSqk9lFI7AC8CfxcRe6XZQ8A4YLxSaizw\nKtAkIiNy5rsZOBaYqJQaB9wFzBGRbulWIjIV+ClwmFJqZ2Aq8KCIHOz6KnMolAbakUq7cAv1YMG4\nsIL8Noks5NKJx2PEYsXt++G/zYu/67DkiISCcajFSSbivLf4M3b57jV5DTmD4MITD2J1cwvNZnGu\nU4zQjcJ3w8F778jsl/O3V9ZpypogKEbBfApMU0otth273nzvUQAiMgk4GLhCKWXl+M4A4hhKAnPc\n9sDpwHVKqVUASqk/AIuBa2zjBgGXAzOVUh+Z454B5gA3uL7KHOLxGNmsylssOlz2zgrCteQ0h986\nmGKz0Iw5SuPeCmKOdCZb3iB/AcshEY/xwJzXObpx55LsrtmvroY5M89lYP96U47yKdvdv7oVf7z6\n5Lzjul2/Jgh6/fYqpR5VSt2Zc3iQ+dtyHE8GOoC5tvelzL8n2953jPm7KWe+Z4GDRMSKiB8K1BcY\nN8ZUVJ4REcfFJBj3lr/UXCg+xbhQ9parWFIimDRlP/vBWHNE1YJJxGPcP+e1orILg5LDHgv6fOVa\nTr/mjyRLoGzj8Ri7fGVLRxm0gtH4xfUnVkRGYri5XgduMQ/vBCxTSuWuuIuBTURkmG1cFljiMC4B\njLGNs47njgPDFeeLQgomCPdWrhWklGLfM27gk89XdZ/D55N/ofe7u47CitJPs0ulVJ9UMLk9wFLp\nDIjQkc7wtR1HlUUO+2dTKcXJV9/D9lts7LvRpRu0gtEEQdEKxgz2v4+hHGLA0Uqp9ebLw4B1Dm9r\nNn8PtY1rUfmObKdxOMzZDIhtnGcKZ+q4CfIXFxwXEb46agT35WyX25MVVLSLzCFN2V2yQqIk1kcm\nkyUWk6LiQIXmiIQFk87Q3pHi1CMnlsQ9VkiOZSvW8uK8D/jVn55m1doNnHzEhLLfC92uX+OXohWM\nUup/SqntMNxj7wPzRWRCySQrMU6LSUfKTe1HYdeSk3I48Zt7cu+T/+4W9ylkBfltFROEqy8I68PN\nE7fTvYiCgkmlM8RjMU45yjHRsSQkEnG+WN3M+b9+kN898C/u+9mpKFO+cqFjMJogcO0iU0qtV0qd\nByyny0W2EhjgMHyg+ftL27h+kv8o6DQOhzlzx3nGWcG4jV0UH4OZMH4b2jtSvL7g424y+Mng6rGO\nxoWi9KPkIJg+Ys7xD+edJEuFVeBofwhIZ7JsMnQgI4YN6uGdwZJMxNlm8+G8eu+lfPzEdWy7xcau\n/k+9sm5DWzcZtILR+KXXT6yI1Cml2hxeehuYLCJJYD6wm4gkcuIwo4HlSilLYcwHjge2oHscZjSQ\nBhbYxgGMchinbK87Mm3atM5/NzY20tjYmDfGqS2IqyyyQgWKBZ7cRYTTvvV1Zj7YxN07/qBzDqfY\nhe8aFpfZcOtb8v973aZbt7Z33/8mCAVTbgsmFosRi0m3FkBGFln5ZIAC8cES34vn/rOIS2c+wkt3\nXgJoBVMNNDU10dTUVNJzFLOCPCUiU5VSr+QcHw00K6VSIvIwcAYwAXgewFQ8E4H7bO95BLgWaATu\ntR3fD5itlLI2N5kFtJrjnreN2x9YoJRa1JPAdgVTCKeAblCV/IUWgh9N/gZHnn8L6bSxv4iTiyuI\nGhb31xF+BlgUFIwlh71FjdN2yWWRoYdWMaVgz7GjeefDZaxu3sCQgQ09tiHSVAa5D9/Tp08P/BzF\nusimi0hn33cR+TGwG3AjgFLqaWA2MENE6s1hl2NYJZ2blpuK4TbgUhEZas51MrA1cJlt3FqMOpqz\nRGS0Oe5AYBJwgfvLzCdZoKGg7yB/DwvBkIENvHDHRZ2bVxWMXRRdw1LIveXyOpwSBdL+4jheFEzY\n7fotOez3tFSNLnuiUD1Osf+nXqirTfL18dvyzKsLgZ47dWs0xVLMY+5PgdMwqvLTGPUpK4HvKqXu\nt42bglGAOc8ctxRoVEp9ljPf2cBVwFwR6cDIFJuklHrbPkgpdb2ItAJPiEgKyABTlFJzXF+lA4WC\n/O4W1ULxD+8V8G6sj8KFmmW+jmSCVNrfE3chZdu/X13RcwSBJYdSChFx7JhdDhkKtesvJYdM2JFZ\nL73LsQfupl1kmkDodSVTSr0MvFzEuBbgnCLGZYArzZ/ext4E3NTbOC84berkqlVMgfYofjO43Cwk\nPdXRlKJdzTsffEp9XQ3bbD6823FHZZ1216iypyLHcmJV0d/4l3+yrqWd/XbfIRQZcl1kbj5XXjlk\n7x25/p7ZKKW0gtEEQlV2U4Yg6mBK04ssmC0Dgumplrug/fO1hXz91F+wzVGXcdnMR1mzzgiZlSqL\nLCwX2fuffME1dz7FlAN2jUyQvxz3YrstN2aXHbZg5Zr1BWNzGo0bqnLDMSj81O2qAj4I95RTLU6R\nSq7QpmX2LLI3Fi7hN/c9w22XfY/6uppu71dKsb6ljXUb2vngky9Y3dzC068s4OPPVjl2dP7JCQfw\n4+P35+0PPuWm+59l+2OuZOpJBzN8cP9eF0QrO662JlnwXrR3pHqcoxwk4jHOveGvXHHaYXx19AiW\nLl8dUqJB4Xb9pUJEeOK3ZwOwurlFWzAa32gFY6PUWWR2Pv1iNU+99A5jt9nM0/u75DCuo5uCSaVp\n70hz7CW3MvetD7n81G867meybMUaLvzt32jrSHHw2csYMrAfO2+/BWcftx8vzPvAUQ4RYdx2m3PH\nFd/n3Q+XcctDTYwcPjjvXn64dAWfrVzDN3/8O5Z8voqPP/+S9o40Zx3byG8uOC5v3qRDunQYCmZ9\nazsjhw/m7OMaO2UoZ8NNcE6hj0LRqUbjFq1gbHSk0vSvry3q/U5ZZJlMFpHiNqcaufEQtttiY359\n3zOcetREhg7u3ymDm4I6S9HV1RqWQTab5dGmefzr9f8y7YdHcPe0H9BQ4JpGbjyEB647nRvv/yez\nfveTbq8Vs6DtuM1mzLzkOzzWNM+xnqeuJsmZU77BViOGstWIodTVJFixer3jXIViMOVcVOe/v5S1\n61q5+kdHdqaJpzPZsm56BoXa9Zc+BpMrg1YwGr9UbQzGudmlCxeZQxdit4vA7mO2YsetR3DYuTez\nodXY5cBtQV2um+3l+f/jP+99zOHf2ImpPzikoHLpvI4C3QDcpTrn90QbPmQAozYbxhHfGM+47TZn\nUP96amuSbL7JEMc5lq1YwxermjtlWbV2Ax9/9iXxWH7/r2w2y4dLVxQlmxt22GoTRm82lGGmsoeQ\nLJiQYjB2CqXhazRuqGoLJu8pMe2vfsRL/62D9x7DgsWfc/SF/8fvLjreMdHgnQ8+paWtg403GsDw\nIQO6KY3cTLKJO2/Led85kDf/+0nRMoRdaJnNZnnujUW8v+QLBnzjxwwfMoC161upr006bgy3eNmX\nTDjlesaMHsEvfjyZPQLqclxbk2Rg//pun4soFVqWsg7GSQZtwWj8UtUKxm8dzLqW9m7HvFgf6UyW\n2y8/katvf4JPV6whmYjnKalHmubxaNM8Vqxex4o160nEY2w6dBC3XHKCYywoqGy4cimYWCzGjybv\ny8vz/8fMS07g0xVr2GrTjTh26m0M6l+fN36bzYfzyZM/58+zXuWI82dywkF7MOPMI4uqmVn08XLm\n/HsBInDWcfv1ei1RySIr1Y6WTrS2dXDzX5u0gtH4RisYG+73g2npdsxNcaIxR7xzIZ5x5lEA/PPV\n9/IWkivNKCK+AAAgAElEQVROO4wrTjsMMGIb6za08fmXzWy80QDHFGE3ii64jszBtIqpq0121tn0\nNEdNMsEPjpjA4fuM46IbH2LH46bzwh0XseWmG3Ubl81mefd/n/Fo0zwefOY/rFyznm9OHMsR+zhv\nKZT7uQiv0DLfgqmvdc7AC5q62iS/f+i5vKw+jcYtVatgCm1L6yeLzG0MJhmP09qen5rbk/UhIgzs\nX9+5vW6hNvfFx5KcW8W4uxel6UVWjHtq2OD+3HXVD3j+jUVsPCS/ofe6lnaOuej3HDphR2ZecgIT\nx2/TY5+3aLSKcb6fAxvK09VARDjxsL06iy7LtQ+OpvKoWgXj9JRo7AdTfBW901O72xhM84buqblu\n9qQBZ+shqJ05o9Dsslj31Dd2dd5Fe1D/et5/ZEbRcuS6p6LUKqacO1qecfQ+XH/PbNa3tDOgTIpN\nU3lUbRZZKepgPDV4dEgU8FIHY8dIdXbTtDM/xdjNwlqokt9NUNqp7U1o3ZRtDx5RaXbppvloEGy9\n+XBiIvzt2TfKdk5N5aEVjA2/WWRuFJQ1hx8lB84bn7lKty6gKBPxWNGukSB2o4xSq5huFkwElByE\ncy9qahL84bG5ZT2nprKoWgXjuw7GIb3XWASKv6WO7i2XT/6O7WZCaLUfhILxa80FgVOQv+xZZIn8\nvYrCuBf9apP85vxjPb1XKcUnn6/SmWhVjo7B2PDrInObShrEwlwowO6uSDII5VCKIH848Y+8IH9E\nCi3LWQcDxmdjs+GDXb/vjYVLOP/XDzL/g6W0tHUwasRQxmw9gtsuO7FbEaum8qlqBeO3DiaIIL9v\nF1kBRecmiywIGUoX5A83gyuMQstEPEY6k+2WwVXOOhgLr8WW9z75Mt85ZA+eOfJc0pksH3zyBctX\nNWvlUoVoBWOjI+0iOF5gYXb95O+Q6uzaReY4R3mTFYK4F1GJweQG+Qt1gC4VsViMWEzIZLKdu5+G\nEoMpUITbG7+94Nud/04k4ozddiRjGRmkaJo+go7B2HBdB+PUKsZlo8pSuMjcpFtbjTkzOYtqGDGY\nqCiYsNOULTlyFV3Y90LHUzRuqVoF49yx1k3swl8FvSVDfg2LSyXl0HTTUwzFdi1uXGyGDNFIFAiC\nZM5Op2G46SDfVVfuOhjI/z+ZcvGt3PzAvwAjiP/6go94c+GSssqk6VtoF5kN964lf4FYRzeby3oH\nvy4y6Ar0W12/3C7s8XiMbFaRzWY7q+RT6Qz9cjY464lCQf5EmRVMfgwmTAumS45yt+u3ZLDfi99e\ncByTzvot/3r9v/xn4RJqknGmnXEEu3xlS9dzz/vvJ3xl1Kad20xoKhOtYGxEpdCy2D1pLDn81PNA\nvqJzex0i0nk/a2u6FEzftGDCbxUD+S7ccGIw3S3s0SOH8cIdF3H/7NeY9sMjGLvNZp7byFw68xGO\nO3A3Tj5yYlDiaiJI1brInPo9daT9tYoJIrDt1j2VTMQKbLvsPRPNbT2PIUf3+9GXFUy6mwVTfhk6\n5ciJwYTtIgMYMWwQ5333QHbadqSvHmXnnnAAv/3Lsyil/IqpiTBVq2CMhcRnHUwA+8H4D/I7pzr7\ni8G4c7E5ydGXFUwqk2vBlF/B5FowfSlNuRgO2msMqXSGpv8sKsn8mmhQ1Qom6DoYbwt78HUwXhSd\nH+sDnF1LvgstI+AiS6d1FlkpEBF+cvz+/PbP/yzJ/JpooBWMiVLKlYLwW0FfaA43KcZQKMjvbV8a\nC68KJt/N5l3BWA03w240GaUssrLHYBz63AXJiYftxUvzPyzJ9teaaFC1CiZ3IclkssTjsc66kN4I\nqlWM30SBQplofhpmdris5wH/Ssqp/iTuouFmUESqDsZ+PzPuHhqCkqGUtS/96mp44LrTGTKgX8nO\noQmXqs4is7sg3AT4oQfXlM9CSy8uMsc9ZVxaUh0BWDBBKpgwntgtOeyxuXQmU/ZUaYBEItb98xmB\nNOVSsP8eXynp/JpwqVoLJvfL475NfjANHvOLNT20inHYNsCPHEEpGDf3Mx6PoVRXR4EwFUxUFF3Y\ncpRDwWgqG61gTNwqGKu4MNPNCip//63CsSDviQJu+6GB6WYLqJYGwgnwQ6F2/eHUwaQzIcdgHDpt\nazRuqFoFk8hTMO7cSiLiGHdw55rK/wK7r4MJIhMtmCB/kHOEFlyPTKFl9FrFaDRuqVoFE0QrDqce\nXv4LLd3J4Zhs4DaelAg/TTl3jqi4pqJSaBlOHUx+Q9hSkc1m2dDaXpZzacpHVSuYbhaMp8ypfLeQ\n2zhOEO6tIBIFcmMw7gstg43jhBrkz6k/iUKhZWgushKmKduZfvsTXHfXU2U5l6Z8VHUWWX4Mxm3c\nIX8Ot00iM5lstyaR3upggt2R0qs157fyPCoKJjJpyqais+J8xabQBylD7sPLmnUtXHTj3wK3bL5Y\ntY6X3/6QpV+s6Xb86MadOapx50DPpSkfvSoYERkPnAV8HUgBceCfwAyl1ErbuAbgeuBAIA0sBc5X\nSi3ImS8BXAVMMedrBi5RSs11OPe5wOnmuLR5zsfcX2Y+Tk+IXtqj+HnyN+I4hgXitUlkrpLzshgF\nbX0EMUeYC3vUYjBuY3tB4aRg3l/yBc+8+h5XnX54oOfKZrPMfesDxozelI03GgjAi/M+4OF/vakV\nTB+mmNXwAeBtYFelVJuIjACeBQ4WkfFKKctx+hDQAIxXSrWLyNVAkznmM9t8NwONwASl1CoRORWY\nIyJ7K6XmW4NEZCpwPvA1pdRHInIg8A8ROUIpNdvndTv4uD1YMI5NIr3Ecbp2TXRdJJnTEy2MPmIQ\nfKpzVFxTobaKMeUIowbGkiFXwbR1pBg5fDA/OGJC4Od7+e3FJOLxzrnrapI80vRm4OfRlI9iHs2y\nwMVKqTYAU1n8EtgO+CaAiEwCDgausCmcGRjWzk+tiURkewyL5Dql1Cpzvj8Ai4FrbOMGAZcDM5VS\nH5njngHmADd4vNZuOLl0vMUdgnMLdcnho1GlywA/5FswbmUw5ghGSXV3kZXfcnC0YMIotLQpurDc\nhU4xmNb2FPW1xe/z44aj9h3PY8+/1fl3fW2S1rZUSc6lKQ/FfIPHKaUW5xxbBggwxPx7MtABdLq5\nlFIp8+/JtvcdY/5uypnvWeAgEbF6RhwK1BcYN8ZUVL5wjMG4DfInHFxkPhd3v3UwniwYpywyl9ZD\n3sLsIfsqKjEYe3ZhFLLIohKPAmht66C+RJuEHbDHV+hfX9tpudXX1dDarhVMX6ZXBaOUckoj2QHD\nsnnO/HsnYJnD2MXAJiIyzDYuC+Tus7oYw103xjbOOp47DmBcb3L3Rl47dA9+7twWK24LLSGAjgIB\ndUKORgwm/Kf2yMRgbPeiI4QaGCigYNpT1NeVRsHU19Xw5I3ndFqM9bVJWts7SnIuTXlw/c0RkRhw\nCnCHUupD8/AwYJ3D8Gbz91DbuBaVv8uQ0zgc5mzGsJyG4hOnGEwQC7Pf+IdbBZNfi+MlBpOvKP26\nC/0H+aOhYKKwZXK4yrb7nkmt7R3U1ZRnm+O6mqS2YPo4Xh7NrsRwh50XsCxlxa/lAM7KwfWiGvfX\naDJfhvDjJ0HMEZkgf2itYnJcZCHci5pEfq+8UsZgcjEsGK1g+jKuvjkicjJGevEhSqlW20srgQEO\nbxlo/v7SNq6f5PdgdxqHw5y54zzjN7gOzvETL3PkF3x6r4PxFEsKKBuuEiv5oyBHFGSwMBRMeSwY\nIwajXWR9maJXIhE5EcNq2U8plbvAzwd2E5FEThxmNLDcVi8zHzge2ILucZjRGHUuC2zjAEY5jFO2\n1x2ZNm1a578bGxtpbGzMG+PY68lTkD/Y+Ie3rs7ds8iCcfV5aHZZgZX86Ux4lpQlh5cU+iBwTFMu\nYQwml/raJG3agikZTU1NNDU1lfQcRa1kIvI94CLgAKXUCvPYYcBmSqnbgYeBM4AJwPPm60lgInCf\nbapHgGsx6mDutR3fD5itlGox/54FtJrjnreN2x9YoJTqcSNvu4IphFMvMr+Fln7TlDOZLEq5LZLM\nd2956UgQdLp1JRVahi1HFGSwaG3vKLmL7Ms167nj0Rc545h9tIushOQ+fE+fPj3wc/S6konId4Hb\ngLsxUom/ax47AhgBoJR6GpgNzBCRevOtl2NYJddZc5mK4TbgUhEZas5/MrA1cJlt3FqMOpqzRGS0\nOe5AYBJwgY/r7cSxF5nPLDK/3QC8KAdHF5kXGdIBWB+mwlbK2MbAbezCbklFZVENtV1/FOpg0k4x\nmNJaMP371XLtXU/R0tahFUwfp5iV6CagFqO4Mhe7ypuC0SpmnohYrWIac6r4Ac7GaBUzV0Q6MDLF\nJiml3rYPUkpdLyKtwBMikgIywBSl1JwiZO4Va4Mrqw9YMBaMvzRlLwuJcx2Mvzm8Kco461raO9+f\nTMRdb3dck0h0uxdht2ixFGW5e4B1ymHLIotUmnKJFUxtTZJ9d92e5998n1Q6E9r/gcY/vX5qlVJF\npQSb7q1zihiXwchEu7KIsTdhKLiSYGUM1dbEPC/MwXQh7moJEoSbztscwTW79PrEHRW3kBX7sKwX\nt4oyKDla2owAd5RaxZQrTfmQCTsy66V3qa1J0NaRoqG+tuTn1ARPVT8W2L9AXuMO3RZmn3EHL21e\nHIskPVThRyFAHxUF082KCkEGiEqrGKdK/hT1daVPUz5k7x2Z/e8F1NUktJusD1P1CqYzUyftrVVM\nfiaa++wrPwuJdQ1W7ap3Kyi4bDivT9xRC/KHtW2zJUcUCi1zt4IoV5ry1psPZ1BDPYl4XGeS9WGq\nXsHYLZggXGR+F2a3yiF3L3tPyQq5qc4e3WwVacGE5PtP5MgRmTTljvLVwfxpxin071era2H6MFWt\nYOxuCO/ZV8HFLrzIYM3hR0mVphOyfwUTxuIei0lncD8sKwpyrOsQtku2ZAgjTdlijx1H0b++VrvI\n+jBVrWByF7RybzjWOYfPhdmuIIJKVoiCggljUbUswnQmE1q7GohKDCZ/O+5yVvID1Nfqjsp9mapX\nMPanRPfxk66F2UuRpCWDHzedNYflK/fupvMZS7LFcfqygrHLEVbDTUuGKLTrz+tF1pairowKpq42\nQWubdpH1Vapewfh2kaX8+cnt1oOXPmL5c3hoeZOb6uxz2wFfCiYTbjflTjnSmdCKLMFqFWP7bHr4\nXPilsItMWzCa4qhqBWMvqjNa1Hup5A8ucyo8F1kwWWTBx2DCdU9FwYqCaCQ8WJSzmzIY/ciaN7Sy\nvqWtbOfUBEdVK5gg62C8Vlt3Uw4erSB7y38vm1PVOPTf8jKHpWwrx0UWpgUTJ50O10XmFIMpZxYZ\nGBbM3Y+/zO8e+FfZzqkJjqpXMOm0PxdZEE/tfutH7AVxQfQi89vsshIUTDqTDV2GsO9FMtG91x5Y\nO1qW14IZu81mPPjMf8p2Tk1wVL2C6R5g9x678FLFnyuDHyvIT5NIpyB/tSuYVDpDOmQ3XVcRcDTq\nYJRStHekqaspXzyovi7JFpsMYdnKtby/ZHnZzqsJBq1g7DEYT1lkNuXgxfpI5PYi828Fuc8AC2Lj\nNH9KDnKr6MOtQdExmHwF09aeoiYZJxYr37JRX1tDeyrN5P130VZMH6SqFUwiEfNlPXSzYEIN8vu3\nxOwLifeCUX8xg9yYVpjxj84YTCK8Sv6wW8Uk4jEymSzZrPH/2tqeKkujSztWL7LjDtydB57+T2dL\nJE3foKoVTPc6GH9ZZL7cW2nvKcbQfWH2lg3XpSiVUtpFZrNgwkxTDvte5LYhKneKMnSlKX99523Z\n7Stbsm6DzibrS1S3gokHEYPpcm/5XVS9torpbkl5jMGku3bVjMXEtRskP5bkzxILW8GkM9nQa3G6\nP/yUvw7GkqNLwZQ3wA9d2ybH4zHuvOokBvav7/1NmshQ3QomkEJL77tRGnPkLqru/0tyOzK7vQ7L\nFZIxM6e8KTn/u1F2i8GE3Aesy4Kp3lYx0N192lbmNjEA9XW60LIvU9UKxh6D8bqXSxDtUbrcW94W\nd7sF4rUjs7WQBJUNVwkusrCVXNitYiw5rM9nufuQgWHB6G7KfZeqVjBWU0PwmH3ls1DTksG/iyy3\nVYx3SyoI5RDEvUilQ9zsKxGzWTDhJhpA+Aqmewym/C4ybcH0XapewfiOwQQR5PfpWureDcBfHMdz\nunXAacpRWFTDt2C8fzaDlMMegyl/FpmzBfPZyrU64N8H0Aqmm/XgpX7EFuT34K/Pk8FnNwBfrrp0\nJtR068gomLjR4SEqhZbhxmBygvwRicH84p7Z7HvGDXy4dEVZ5dG4o6oVTK4bwn+Q36vl4O9JNb9V\njNeCz3SoyqF7kD98CyYKVhSEb81Zn/HWtg7q60KIwbTlK5hfn38sJx8xgb1Pvp6//fONssqkKZ5w\n7O6I4DcVtBStYjy5yBJBKCnjWtKZrOdEg3Qm67mOxpojKotqZ6FliFsmd/bJ89BlIijy0pRDiMG0\ndeQrGBHhnOP3Z6+dtua4S2/j9kdf4OQjJvDtg/bIG/tY0zyee2MRozYbxl5jR/O1saPLIboGrWBI\npTNks1lPi0mNLXvLc5pywtby30ORJPjfMhm64jheF3YR8d3mPi/IH9rirtv1O8lR7k7K0Pt+MHvs\nOIo377ucZ19byOjNhjmO2WLTjdhs+GAWfbycX9/3DOO2Hcl1Zx/NjttsViqxNSZawZgLSU0ygYi4\nen9kW8V4nKMj5W8XR/v99LIQRW1RTWcy4bWKiUoMxlZjFdU05cED+nHM/rsWfH3Xr2zJrl/ZEoBf\nnTeFWx58jv3P/DWv3D2VUQWUkiYYqlrBWE/cftJqg9wm2Kt7K5ljSXm/ljSZbNZz5157JtrAhjpP\nMkRFweh2/V1yWJ+t1rYQ0pQLBPm99iSrSSY49zsH8MNj9qG+rqboedw+fGoMqlrBWAuJny7GqXSm\nM+7gP8jv3b21dn2rMYfXNOVEgg7TXRiEBVMJlfyG2zT8LLJIpSmX2YKpq0nQ2pZvwZx45Z3c99Sr\n5ZGhNsnns3/JIN2mxjVVr2Ba21OeK+jj8RixmJAxlVR4LrJEjhXk1UWWJptVkVAwYT+16xhMvhyt\n7R0MGdivrOcvFIP5+LNVPHfbBXxj1+1LLsO237qcFavXaQXjgapOU7a+PF7rT6Ar0O+nF5nfKvwo\nBPktOfwqGHvad9hV9GFmkVnNRrPZcF11YdfB1NYkSJt98uysXd/KwIZgF/wF/1vm6DIb2FDX6SHQ\nuKOqFYw968mrC8KKOwTSKsazeyun4NPH4m5U8vu5F5VlwYSlYADfWXlB0K0OJoQ0ZRGhriaRl6rc\nvKEtUItCKcWpM+7lt3/+Z95rg/rX06y7BniiqhWM1Y7DTzv07i1Wwm8V4z0WZCgprwoKKisGk874\ny6gLTg7vMcKgZOhMUw7BggHnQP/a9a2BKhgR4S/XnMbP75nFC2++3+21Qf3rtQXjkapXMFYWmffM\nKXNhTge1H4zXNi9dBZ/etm7uUpRRUDBhP7V3WTDhKZgoWDD2h5cw9oMBq5q/K9CvlKJ5QysDPGQq\n9sSozYZxz7QfcPxP7+CzlWs7jw9s0ArGK0UrGBEZISKzRCTb++i+QfcYjH8LJqw0ZctPbmxvq4h7\ncOsYjTv9ugttbraQugEEQWcWWTo8K8qSI0rp0mHsaAn5gf71Le3U1SRLck8OmTCW07/1dY6belun\nZ2FQfx2D8UpRK5GIHAO8BGwNFEwcF5EGEblZRBaKyDumQhrjMC4hIjNE5D0RmS8iL4rIxAJznisi\n74rIPBF5XUSOKvLaesUK5npt8wJGY0R/Qf6E7xoWe+yjJhn3lLNvKSlf98KnBROLxYjHY51fbLe7\nagaFZTmEGeSHrnYxYSsYewym3N2UwUxVthVbNm8I1j2Wy5WnH8bWI4ex8KPPATMGoxWMJ4r99lwE\nHAjM7WXcQ8A4YLxSaizwKtAkIiNyxt0MHAtMVEqNA+4C5ojIOPsgEZkK/BQ4TCm1MzAVeFBEDi5S\n7h5JJuKk0/583J0Ls8cgv7WbZDab9RwLsgfo/SiHoFxkXt2F1hwtbR2hWw5hu6bsckRny+RoWDBB\nx19yicVi3DP9ZMZttzlgxmA2aAXjhWIVzESl1Ic9DRCRScDBwBVKqXbz8AwgjqEkrHHbA6cD1yml\nVgEopf4ALAausY0bBFwOzFRKfWSOewaYA9xQpNw9EsQXuFuQ38McIuJ7QevcjTKw6/Cast3lZuvr\nCiadyZqtYsKPwYSZ8JCXphxGDKYuSVuOgvHSKcIrOgbjnaJWI6VUMXGXyUAHNitHKZUSkbnma+eY\nh48xfzflvP9Z4Ici0k8p1QIcCtQXGPdLEdleKbWoGPkLEVyQPx1I9pXXAL3Vssbr+6ErWSEeD6/Q\n0pqjtT0VuoKJigXT2p4iHo+F1qrE3g4ptCyynF0tg05R7g0vWWTZbJZ5i5byz1ff45lXF/LOh8vY\naGA/Zt/8EzYbPrhEkkaPIO3unYBlSql0zvHFwGEiMkwptdIclwWWOIxLAGOA181x1vHccWC44nwp\nmEQ8RiqTCaz2w5/14N0CsToy+7JgzCyyeDYWuoKJggUTiRhMPDr3AsIptATLRdYVgym1iywXqw7m\n1XcWs8eOo4pS9kdf+HsWfvQ5k/b8KmdO2ZdddtiCtetbGT5kQBkkjg5BKphhwDqH483m76HASnNc\ni8ovmbWPs+bDYc5mQGzjPNM9BuNvYfYag7HkCN9FFieVyZBVyt8cASmYai9wBMuC6QhtLxgwklgs\n6yGMQkvIt2BKUcXfEwMb6li9roWzf3E/o0cO5Zc/mcIWmwzpUdHcPe0khgxsKGr+bDYbWkJLqanM\nqyqSbq4pPy6yABZVS0l53w8mCBmCuY6Uj/TeaFkw4Vfyh30vutfBdJS92SVAXY2Ti6x8MZhB/etZ\nt6GN52+/kOGDB7DnSdcxqPFc9jzpOm78S37lP1C0cgE47NybeeaV94ISN1IEacGsBHKzxQAGmr+/\ntI3rJyKSY8U4jQMYAKzuYVwe06ZN6/x3Y2MjjY2NjuPsQf4gsq/8xD/8WCBdcSA/ijJBR2oD2WyM\nhvpaT3ME4iKLiFso7PoTS44o3IuOdJpsNkt7R5q6mvJns9XXdS+0DMtFVleb5OZLTuDmS05gdfMG\n3lv8OZms/7LAqScdwrFTb+OF2y9kh1GbBiBxcTQ1NdHU1FTScwT5aZkP7CYiiZw4zGhguRl/scYd\nD2xB9zjMaCANLLCNAxjlME7ZXs/DrmB6IhFgFpnf1NzAXGQ+C0aVioceg4lKkD/Mdv3QZcGElaIM\nXfeirSNNbU0iFFeOU5ryqBG+PeRFM9AhyD9kYAMTxm8TyPz77rY91531LY44fyb/vmsqGw0q3vrx\nQ+7D9/Tp0wM/R5CfloeBJDDBOiAiSWAiRn2MxSPm78ac9+8HzDYzyABmAa0O4/YHFvjNIIPgssh8\nN8xMJGj32Q3AUHJ+UozjpNJZn5ZYMO1mNrS2hx5cj04MJlxla32+w9hszKK+Ntmt2WW5LZgB/WrZ\n0Nqe19E5SE791tc58hvjOf6nt3veTC2KuP0WF4xqKaWeBmYDM0TE+t+/HMMquc42bhFwG3CpiAwF\nEJGTMboEXGYbtxajjuYsERltjjsQmARc4FJuR6ymhr6e/BNdGWBBxB28VeEHVwcTZrKCNUcU3EJW\nq5iwFV1U7kVbRzgZZOCUplxeBROLxehfX8u6ltJ2VL7+nGNYvmodz/3H97NzZChqNRKRX2As7FuY\nf79hvvS1HHfYFOB6YJ6IpIGlQKNS6rOcKc8GrgLmikgHRqbYJKXU2/ZBSqnrRaQVeEJEUkAGmKKU\nmuPmIgsRTJA/0RnY9mMFbWht9xcHCiLRwGyYqRWM8eCRSodcaJkIP8hvJX+ElaIMhovsi9VdyaTl\ntmCgKw4zeEDpNlyLx2M03Xp+Sc9RboottLy4yHEtdBVU9jQuA1xp/vQ29ibgpmLO75Zg9oMxFne/\nFsyGVu++9ppEV7Gn367QgvhsdlkZCsbKIgtVDtOCCatVP3Tdi7BSlMHqphxeJT+Ur5rfTfZZX0Cn\nKQfYYsVv3MGPFdUZSwqpKzRUapA/zGaXsdDvhfXZam3roK42nGSDupAr+UHvCeOVqlcwnRs6+Wqx\n4r0XmTFHgg0+nto7W8X4VJS+G1UGEBzvKrQMuweYTlO2ZIiEBRNiJT9oBeOVqlcwfjohW3MEEeQ3\nLBhvysHoVYWvp117oaWveFTanxVUk4zOohq6BRORIH9HKh1aJ2UwYjC5WWQDy6xgBjbU6Zb9Hqhq\nBWPFYDrSPlvFpP1aMGaQ38dTe00y4c/NZgb5/Vofbe0pRPC06Zk1RxQWVSvIH7Ycre3RqIMJ34Ix\nFEx7R4pMJlt2ZReGBfP4829x6c2P9D4wwlS1ggmyDsbvRl1GkN/7YpZMxFnvwwoKKk25pc1fzMBa\nVMNe2LtaxYRfaOnnwcMvlus0rE7KAPV1NZ1Bfiv+Uu7u0mEomN3HjOLWh59n5Zr1ZT1vkFS9gumM\nwfjshOyrQDGRYEObd+VgzBFnfYufRAG7i8ynJeZTwUTBgolKoWVU7oWxF0yYdTBGDCaM+At0pSmX\nkxHDBjF5/12Z+dd/lfW8QVL1CqbLgvH75O8/TdnPQmIkCrT7yiKLQg1LV5A/3K2KoxGDiVIdTEco\n2yVD92aX5e6kbDGwoS6UIP+F35vEzAefo8XWi60vUdUKxuqr1NaR8mF9xGnrSJHNKl9xBz/WhzWH\nHyXV2W7Gb0cCn+6tKD21R8GCaW1PRSgGE34lfxgpyhBeFtkOozZl4vhtuPOx3narjyZVrWDAfwaX\nEVw3ArFe/cKWa8mXi8xvkD8ZzHbHwbjIwq+DSWeypMOu5I9AFllNFIL8dTWRcJGFlaZ8yUkH879P\nV4Rybr+E92gUEawnZj+xi/Ut/hdVvzGYLivIXzeAmkTcl7KtqCB/BCyYsBVMMmFsRNfa3kG/MJtd\ndrdD410AAA5pSURBVHORlbeKH4yOyuWOwVjstdPW7LXT1qGc2y/aggnAtWQoB3/xk5a2Ds/xE2OO\nILLIohODCX1RjUIMprMXWXgyWK7TtvYU9XXhpyk3V6EF05epegWTiMd8WQ+Wa8q3BeMzTbkmYbrI\nfBQ4BpNuHcy9SIS4qCbiMTKZrFnAWt50WDvWdsVhxmA6W8WEGIOpSSZIZ7JkMtlwXWQbtIJxS9Ur\nGN+NJm0xGK/UJBOsD2hh9hvk97sfTCVYMCJCPB4L1XoBI5utLSJ92cJUMCLSmaq8tsyt+i20BeMN\nrWCsGIyPLLJglIPfIH8whZZRcJEppUJdVC05wiyyBDoVXDQUTHhpytDV8LJ5fVs4MRizVUwUNgMr\n5cZnQaMVTGeQ36eLzMdiFFgdjI8sMitzqr3D/9bPfhUMEGr1uiVHFJSc/XdYMnSkMrS2hReDga44\nTFgusppkgoSZNh4mz762kKMuuCVUGdygFYz5xOxvszC/8RP/+374bRUjIv6VbQALYhQWVev8obvI\nzPNXex0MmA0v21OhucggGm6yieO3Yd6iT3hj4ZJQ5SiWqk9T9vsltlxLQbiF/GaRtfkMCBsBXX87\na9p/+5kjzPoTMD4X8Vi4CiYKytbq1L2+tS1kBWPGYELopGxhuclGDBsEwJsLl7D7968lmy2/22y3\n711T9nN6oeoVjPXl9aNg/Lzf/l5fLrKE/zmSiXhngNvr+4OQwe8cQZBMxImFmEEGdMaAonAvmje0\nhVZoCV0usub14VTyQ74FM2/RJ3zn4K/xxxmnlFWODa3t7DD5Su6/9nS+vvO2gc0rcltgc1loF5nP\nBS3IRdWviwzwnKxgnd9vPQ/4i59EScFEQQbw938alBzrNoRtwdQYMZgIucj++/Fydthqk7LL0VBf\ny6/OncL/+/mfSaczZT+/G7QF4/NLXBOAcrDeG4QV5HeOYBRluNZcECQTcYSwLZjws8jA+D9p3tAW\nbpC/LklrW0dolfyQ3/By0ZLlnHDw10KR5bhJu/PZyrW0daToH/LnoyeqXsFYbgi/LrKw3UJBWEE1\nPp/aK81FFraCidK9WN3cQl1NeMtFXU2SDW0dbGhtZ0C/cBRMbsv+/368nO233DgUWUSEc79zYCjn\ndkPVK5hkIk4s5j3uYC3o4T/5ByNHFGQAIpDBFSfkEExnokMUFIxSKvQYzIrV62ior/X8XfWL3UWW\nyWT536cr2W7L8rvI+hJawfheVE3XlK8MMGsO//GPKLjIwp4jCJKJcNvEGDL4f/AIUo6wYzCff9kc\nWvwFuiuYJZ+vYvjg/vQL0W3YF9BB/kQ8kIU97EU1KCsoKtcRvoKJTh1M2PfC+n7UhZymvPzL5tDi\nLwADG7oUzH8//pwdtto0NFmciGKFf9UrmEQ85mtRDqIYLhjrI4gsskRAGWDeP1ZRUjBRkMH+O0w5\n6mqToVp09XXJSFgwVgxm0ZIv2H6rcOIvTiil2PeMG5j98rthi9KNqlcwfl1kIhJg9pX/OhhfSiqR\n8CVDLBYjHo9VkAUTfrGnJUuYJBPxUN1jYLjIlq8KX8FE1YIREX7x48mceOVdPP78W2GL04lWMAE8\nqRpKyl/2ljWPHxmCmCOIe1EpQf4oyAD+HjyCoCaZCLXRJUBdTYLPQ3eRdaUpL/r4i9AyyAoxYfw2\n/O0XP+SC3z7Efj/8Fc+/sShskbSC8ascILjYRSAuMp/1OFFRMFF4ao+CDPbfYcoRBQvm8y/Xhm7B\nNG+IpgVjsc8u27Hgr9P4wRETOHn6PbzzwaehylP1WWSJRMxXBhhEK/vKbyzIb8ZSFGppgiCZiIfS\nY8pOtFxk4WZL1dcmae9Ih65g1q5vpaWtgxVr1rPlphuFJktPJBJxTjp8b7536J4FU7pHHnoJmWwW\nwXBtX/Dd0tTUVL2CCcSCCSrVOfQ5tAVjlyNsBROle1FfF7IFY6YDD2oIX8G8v2Q524wcHlo9TrH0\nJN+7f72K1vYUSimUUvTvV8cFJZBBKxifygECbLESQLq0n7hBTcLfdUCQMZjwA+xZCdmCiUgdTE0y\nCi4y4/xhdVKGrhhM1DLIvDB4QD8GDyj9eaKtggERGS4ifxKRhSLynog8KCIjg5o/mBhMwneTSUsW\nr1gLu59U0qDuRaVYMFGQwf47TDmi4CIDQnWRNdTX0pHO8O6HyyIZf4kikVYwIpIEngGSwFeBMcAG\n4F8i0i+IcyTiQVgwEQjyR8ASg8pykSV81PMEQbRiMGFbMKaLLEQFIyIMbKjjtQUfRS6DLKpEWsEA\nPwDGAhcrE+ASYGvgzCBOYFTy+w1s+wuOd8VP/LnIws6Gg8pSMFGQwf47TDlCT1O2XGQhpikb56/n\ntQUfawumSKKuYI4BliilPrYOKKWWAwuAyUGcIEq1H37nCMKC8Z9RF/69CIIoFVqGHoNJJMIP8kfA\nRWadf8XqdWwfwj4wfZGoK5hxwGKH44uBnYI4QSXVweTK0NTU5FqOYIpO/bn6IPhCS7f3IuGzI0EQ\nlErZevlcVGoMxu29GNS/jo0GNTBscP9A5ahUoq5ghgHrHI43A/1EpNbvCfz2IgP/7ikR8S2HU6KB\n2y9PJbvIvCyqUankj8K9qNQYjHsFUx/KLpZ9FZ2mHAELxpIjdBdZBNKU4xEKbEelXX8U7kX4CiY6\nMZihg7T1UixRVzArAads7YFAi1Kq3e8J6muTvs3/fnU1vudoqK/1tWOgIYO/RaBfXY1vX7tfOUSE\n/v1qQ4871NfWkAp5v/NEPEZtTSJ0S6pfXQ39Q9pF0qKhvpb+/WqpDTnZYKOB/Ri58ZBQZehLiJGY\nFU1E5ClgB6XU1jnH5wPrlVITHN4T3QvSaDSaCKOUCtRsj7oF8zDwexHZUim1BEBENsGoibnE6Q1B\n3yCNRqPReCPqFkwSeA14D/geoIA/ABOAXZRSLSGKp9FoNJoeiHQWmVIqBUwCMhi1L+8C/YH9tXLR\naDSaaBNpC0aj8YKIvABMBEZZrlWNppoQkRHAXcBBSqnQDIlIWzDFUuqGmFFERMaLyG0iskBE3hKR\nd0TkRhEZljOuQURuNu/NOyIyS0TGhCV3qRGRyRjKJe/JqZruhYhMFpHnROQ1EfnQ/P1d2+tVcS9E\nZHcR+YeIvGt+T14RkSk5YyrqXojIMcBLGC21CloQxV63iCREZIa5ts4XkRdFZGJRwlj7AfTVH4xG\nmG8BDwBi/twNLAL6hS1fCa97IfAgUGf+PQIjVrUQqLWNewp43joGXA18AYwI+xpK9FlYBDyO4Vbd\nMuf1qrgXwHnA69Z1AXHgT8Cd1XQvgK2ANeZ6YHlrfghkgcMq9V4ALwPbYFgwmR7GFXXdwO/NdWUj\n8+9TMZoOj+tVlrBvRgA383RzMdnKdmwTIA1cELZ8JbzuBcDonGOnmPfiaPPvSeaXaV/bmCTwJfC7\nsK+hBPfkPHMhvSpXwVTLvQBGAW3ArjnHN7WOVdG9ONP8HIzLOb4GuK9S7wUQM38XVDDFXjewvXkP\nT8p5/zvA473JUgkuspI3xIwo45RSuX3almFYcFYl2GSgA5hrDVBG4sRcKuzeiMhGwIXApQWGVMu9\nOBFYrZR6w35QKfW57Vi13Iu0+Tu3OjNGV3ig4u6FUipbxLBir/sY83dTzvufBQ7qbduUSlAwJW+I\nGUWUUmmHwztgPJU8Z/69E7DMYexiYJPceE0f50rgXqXUJwVer5Z7sTfwkYgcIyLPmzG6uSJysm1M\ntdyL+zFcO5eb8QYRkcuAGgy3D1TPvcil2OveCWNNyU2WWYxRR9ljrKoSFEzJG2L2BUQkhuEiu0Mp\n9aF5uKd7AzC0HLKVGhHZDpgCXNvDsKq4F8AWGHsoXQBMVkqNAX4D3CYilnVXFfdCKbUOOBCox2g7\ntRxjj6lJSinrIawq7oUDxV73MIy2XLnJAkXdn0pQMBqDKzFM3vPCFiQEfg5cZy4o1U4d0A+4UCm1\nAkAp9RDwGPBTEQl3Q5UyIiLbA69iPG0PVkptDFwOPCIiB4cqXJVQCQqm5A0xo47p/pgCHKKUarW9\n1NO9ASOg16cRkX0wnth/bz/sMLTi74WJpWTfyjn+Jobi+SrVcy9+BgwCzrXWAaXUAxiZU/eYVn+1\n3Itcir3ulRieoNzvVFH3pxIUzHyMzJlcRgNvl1eU8iMiJ2JYLfsppXL/s+cDm4lIbs+50cBypdTK\ncshYYg7E+By/JiJviMibGKmoAP8wjx1CddwLMGIOkP/dztiOV8u9GAssdXjIXAQMx7jearkXuRR7\n3fMxPjNbOIxLYyRTFaQSFMzDwFYisqV1wNYQ86HQpCoDIvI94CLgAMsdIiKHicjp5pCHMTJoJtje\nk8QoRKyIe6OUukoptZ1SalfzZxe6rJlDzWOzqIJ7YfK4+XtczvGdgFaMdkvVci++AEaYloqdURgF\niKupnnuRS7HX/Yj5uzHn/fsBs1VvLbvCztkOIOc7CcwD/oJRUBbDyP/+L5VdaPldoAU43/y39fN7\n4ErbuH9gZJXVm39Pxwh29skisiLvzTRyaqOq5V6Yn/9/Y6SRNpjH9sGojZlaZfdisvk5+Jnt2H5A\nO/CnSr8XGAWmPRVaFnXdwP9hFHEPNf8+GaPQcqfeZKiIXmQiMhwjU2YPjJS6dzD8rp+GKlgJEZEv\ngcEFXp6ulLraHNcPuB44CMOkXYpxb94ri6BlREQOxcgk28T8eQ/oUErtar5eFfdCRAZjXOfBGFZL\nO3CTUupO25hquReTgKkYhaYZjPXhXoxiwpQ5pqLuhYj8AqOQcguMmjgrHvc1ZUtLLva6RSSOUbx8\nHEYi0TrgIqXUS73KUgkKRqPRaDTRoxJiMBqNRqOJIFrBaDQajaYkaAWj0Wg0mpKgFYxGo9FoSoJW\nMBqNRqMpCVrBaDQajaYkaAWj0Wg0mpKgFYxGo9FoSoJWMBqNRqMpCVrBaDQajaYk/H/x+CUzDscc\n4wAAAABJRU5ErkJggg==\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "pyplot.plot(z, p_2016a, color='#003366', ls='-', lw=1)\n", + "pyplot.plot(x, p_2016, color='#003366', ls='--', lw=1)\n", + "pyplot.ylim(0,5000);" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "From the figure we can find that this model do well at the age older than 50 and has error which cannot be ignored at age younger than 50. Because, we choose death rate and birth rate as a constant to make the calculatoin easier, and it is hard to get the actual death rate as a functoin of age. However, according to the data, we can find the death rate in young person is lower than that in older person. That will cause the error. Same things alos happen at the old age. The actual population density at this part is lower since the actual death rate is higher the the value we assumed.\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Reference\n", + "\n", + "Misra, Jagadis Chandra, ed. Biomathematics: Modelling and simulation. World Scientific, 2006.[Read on Google books, page 437-438](https://books.google.com/books?id=U05hDQAAQBAJ&pg=PA437&lpg=PA437&dq=density+function+of+population+age&source=bl&ots=3RkkMcyS_V&sig=jYsJcOvhAsEQKnxQKVF4WG05JRA&hl=zh-CN&sa=X&ved=0ahUKEwj9sKiqyPHQAhUTlpAKHVfjBXEQ6AEIXDAI#v=onepage&q=density%20function%20of%20population%20age&f=false)\n", + "\n", + "Keyfitz, Barbara Lee, and Nathan Keyfitz. \"The McKendrick partial differential equation and its uses in epidemiology and population study.\" Mathematical and Computer Modelling 26.6 (1997): 1-9.\n", + "\n", + "Death rate,( 1 ) United Nations Population Division. World Population Prospects, ( 2 ) Census reports and other statistical publications from national statistical offices, ( 3 ) Eurostat: Demographic Statistics, ( 4 ) United Nations Statistical Division. Population and Vital Statistics Reprot ( various years ), ( 5 ) U.S. Census Bureau: International Database, and ( 6 ) Secretariat of the Pacific Community: Statistics and Demography Programme.[Here](http://data.worldbank.org/indicator/SP.DYN.CDRT.IN?end=2014&locations=US&name_desc=true&start=2006&view=chart)\n", + "\n", + "Age and Sex Composition in the United States: 2006,United States Census Bureau.[Here](https://www.census.gov/population/age/data/2006comp.html)\n", + "\n", + "Birth rate in the United States from 1990 to 2014[Here](https://www.statista.com/statistics/195943/birth-rate-in-the-united-states-since-1990/)\n" + ] + } + ], + "metadata": { + "anaconda-cloud": {}, + "kernelspec": { + "display_name": "Python [default]", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +} diff --git a/population2006b.png b/population2006b.png new file mode 100644 index 0000000..437ab09 Binary files /dev/null and b/population2006b.png differ