{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Foundations of Computational Economics #25\n",
"\n",
"by Fedor Iskhakov, ANU\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"## Newton-Raphson method with bounds\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"\n",
"\n",
"[https://youtu.be/1fD-BSiO1f4](https://youtu.be/1fD-BSiO1f4)\n",
"\n",
"Description: Robust implementation of Newton method for problems with strict bounds."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"- example of how elements of different methods can be combined\n",
" in a poly-algorithm \n",
"- deeper merge than a multi-stage poly-algorithm we mentioned before, such as\n",
" 1. Robust slow method for initial search of parameter space\n",
" 2. Fast accurate method invoked from a very good starting values "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Example\n",
"\n",
"Solve the following equation\n",
"\n",
"$$\n",
"f(x) = a \\log(x) + b \\log(1-x) + c = 0, \\; ab<0\n",
"$$\n",
"\n",
"This equation arises in the models of discrete choice in IO, for example when computing a mixed\n",
"strategy equilibrium in a two players game with binary actions."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Theoretical properties and numerical difficulties\n",
"\n",
"- $ x \\in (0,1) $ are strict limits, any algorithm will break down if stepping outside \n",
"- there is exactly one solution for any values of parameters $ a,b,c $ (where $ a $ and $ b $ have opposite signs, without loss of generality assume $ a>0 $ and $ b<0 $)\n",
" - $ f(x) $ is continuous in its domain $ x \\in (0,1) $\n",
" - $ f'(x) \\ge 0 $ everywhere in the domain, so the function is monotonically increasing, and vise versa\n",
" - when $ x \\rightarrow 0 $ from the right $ f(x) \\rightarrow -\\infty $, and when $ x \\rightarrow 1 $ from the left $ f(x) \\rightarrow \\infty $, and vise versa \n",
"- but depending on the value of $ c $ the solution may be arbitrary close to 0 or 1! "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"- Newton — can easily overshoot to the outside of domain \n",
"- bisections — may take forever to converge \n",
"- successive approximations — ? "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Newton-bisection hybrid algorithm\n",
"\n",
"1. Initialize with interval $ (a,b) $ such that the function\n",
" has opposite signs at the ends, and starting value $ x_0 $ \n",
"1. Compute the Newton step and check if it remains within the bounds \n",
"1. If yes, make Newton step and continue to next iteration \n",
"1. If not, make a bisection step and continue to next iteration \n",
"1. Iterate until convergence, or until maximum number of iterations is exceeded "
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"hide-output": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"{'divide': 'warn', 'over': 'warn', 'under': 'ignore', 'invalid': 'warn'}"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"plt.rcParams['figure.figsize'] = [12, 8]\n",
"np.seterr(all=None, divide='ignore', over=None, under=None, invalid='ignore') # turn off log(0) warning"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"hide-output": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"source": [
"# code up the function\n",
"# make a plot of the function\n",
"# code up the solver"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"hide-output": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"source": [
"#\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"# solution below\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"#"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"hide-output": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [],
"source": [
"def newton_bounds(fun,grad,x0=None,bounds=(0,1),\n",
" tol=1e-6,maxiter=100,callback=None):\n",
" '''Polyalgorithm that combines bisections and Newton-Raphson\n",
" to solve fun(x)=0 within given lower and upper bounds.\n",
" Callback function is invoked at each iteration if given.\n",
" '''\n",
" a,b = bounds\n",
" sa,sb = np.sign(fun(a)),np.sign(fun(b)) # a and b signs, never change\n",
" if sa*sb > 0:\n",
" raise(ValueError('Function has the same signs at the initial lower and upper limits'))\n",
" x0 = (a+b)/2 if x0==None else x0 # midpoint is default x0\n",
" for i in range(maxiter):\n",
" f0 = fun(x0)\n",
" newt = x0 - f0/grad(x0) # Newton step\n",
" if not a < newt < b:\n",
" a,b = (x0,b) if np.sign(f0)*sa > 0 else (a,x0) # update limits\n",
" x1 = (a+b)/2 # bisections step\n",
" step_type = 'bisection'\n",
" else:\n",
" x1 = newt\n",
" step_type = 'newton'\n",
" err = np.abs(x0-x1) # save error for both steps\n",
" if callback:\n",
" callback(iter=i,type=step_type,err=err,x0=x0,x1=x1,bounds=(a,b))\n",
" if err < tol:\n",
" break\n",
" x0 = x1\n",
" else:\n",
" raise(RuntimeError('Failed to converge in %d iterations'%maxiter))\n",
" return x1"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"hide-output": false,
"scrolled": false,
"slideshow": {
"slide_type": "slide"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAsEAAAHSCAYAAAANGxbcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deZiU1Z328e8BGhARBQFxQAQccAGRGESMGo1iVBLjHs1EYjIqMY7rmLwzyZhMlpmJWRyzGaNGY+Ikxj0uo7gjxh2QHY2KCgIKyKLsNJz3j9M93SBLNV1dT1U938911fV0V5Xdd/JEc3v41TkhxogkSZKUJ62yDiBJkiSVmiVYkiRJuWMJliRJUu5YgiVJkpQ7lmBJkiTljiVYkiRJudMmi1/atWvX2KdPnyx+tSRJknJkwoQJi2KM3TZ9PpMS3KdPH8aPH5/Fr5YkSVKOhBDe3tzzjkNIkiQpdyzBkiRJyh1LsCRJknLHEixJkqTcsQRLkiQpdyzBkiRJyh1LsCRJknLHEixJkqTcsQRLkiQpdyzBkiRJyh1LsCRJknLHEixJkqTcsQRLkiQpdyzBkiRJyh1LsCRJknLHEixJkqTcsQRLkiSpZaxbB0uWQG1t1kk+whIsSZKklvHss9ClC4wbl3WSj7AES5IkqWWsW5euNTXZ5tgMS7AkSZJahiVYkiRJuVNfgtu0yTbHZliCJUmS1DLqPxDnSrAkSZJyw3EISZIk5Y4lWJIkSbljCZYkSVLuWIIlSZKUO5ZgSZIk5Y4lWJIkSbljCZYkSVLuWIIlSZKUO5ZgSZIk5c66ddCqVXqUmfJLJEmSpOqwbl1ZrgKDJViSJEktZd06aNMm6xSbZQmWJElSy6itdSVYkiRJOeM4hCRJknLHEixJkqTcsQRLkiQpdyzBkiRJyh1LsCRJknLHEixJkqTcsQRLkiQpdyzBkiRJyh1LsCRJknLHEixJkqTcsQRLkiQpd9atgzZtsk6xWZZgSZIktYzaWleCJUmSlDOOQ0iSJCl31q6t/BIcQtgjhPBkCGFmCGF6COGSuue/G0KYG0KYVPcY2XJxJUmSVDFWr4Yddsg6xWY1ZVK5Frg8xjgxhLATMCGE8Gjda1fHGH9a/HiSJEmqWKtXQ/v2WafYrIJLcIxxPjC/7usPQwgzgZ4tFUySJEkVbtWqsl0J3q6Z4BBCH+BjwAt1T10YQpgSQrgphNB5C3/N6BDC+BDC+IULF25XWEmSJFWIGMt6JbjJJTiE0BG4C7g0xvgBcC2wFzCEtFJ81eb+uhjj9THGoTHGod26dWtGZEmSJJW9detSEa6GleAQQg2pAP8xxng3QIzxvRjj+hjjBuAGYFjxY0qSJKmirFqVrpW+EhxCCMCNwMwY4383en73Rm87GZhWvHiSJEmqSKtXp2uZluCm7A5xKDAKmBpCmFT33LeAL4QQhgAReAv4alETSpIkqfLUrwSX6ThEU3aH+CsQNvPSg8WLI0mSpKpQ5ivBnhgnSZKk4qsvwWW6EmwJliRJUvFVywfjJEmSpII5DiFJkqTcKfMPxlmCJUmSVHyuBEuSJCl3LMGSJEnKHcchJEmSlDuuBEuSJCl3XAmWJElS7qxcma6WYEmSJOXGihXQrh20aZN1ks2yBEuSJKn4li+Hjh2zTrFFlmBJkiQVnyVYkiRJuWMJliRJUu5YgiVJkpQ7lmBJkiTlzooVsOOOWafYIkuwJEmSis+VYEmSJOWOJViSJEm5YwmWJElSrmzYkGaCLcGSJEnKjZUr09UPxkmSJCk3PvwwXXfaKdscW2EJliRJUnEtW5auO++cbY6tsARLkiSpuCzBkiRJyp36ErzLLtnm2ApLsCRJkopr6dJ0dSVYkiRJueE4hCRJknLHEixJkqTcWbYMQvCwDEmSJOXIsmVpFbhV+VbN8k0mSZKkylRfgsuYJViSJEnFtXRpWW+PBpZgSZIkFduiRdC1a9YptsoSLEmSpOKyBEuSJCl3Fi60BEuSJClHamthyRJLsCRJknJk8eJ07dYt2xzbYAmWJElS8SxalK6uBEuSJCk3LMGSJEnKnYUL09VxCEmSJOWGK8GSJEnKnfoSvOuu2ebYBkuwJEmSimfRIthpJ2jXLuskW2UJliRJUvEsXFj288BgCZYkSVIxVcCRyWAJliRJUjFZgiVJkpQ7CxdWVwkOIewRQngyhDAzhDA9hHBJ3fNdQgiPhhBeq7t2brm4kiRJKlsxwrvvQo8eWSfZpqasBNcCl8cY9wWGA/8UQtgP+Ffg8Rhjf+Dxuu8lSZKUN4sWwdq10LNn1km2qeASHGOcH2OcWPf1h8BMoCdwIvD7urf9Hjip2CElSZJUAebOTddqKsGNhRD6AB8DXgB2izHOh1SUge7FCidJkqQKUs0lOITQEbgLuDTG+EET/rrRIYTxIYTxC+vPlJYkSVL1qNYSHEKoIRXgP8YY7657+r0Qwu51r+8OLNjcXxtjvD7GODTGOLRbBWygLEmSpCaaOxdCqK4PxoUQAnAjMDPG+N+NXroPOLvu67OBe4sXT5IkSRVj7lzYbTeoqck6yTa1acJ7DwVGAVNDCJPqnvsWcCVwewjhHGA2cHpxI0qSJKkizJ1bEaMQ0IQSHGP8KxC28PLRxYkjSZKkijV3LvTtm3WKgnhinCRJkoqjglaCLcGSJElqvlWrYPFiS7AkSZJypH57tF69ss1RIEuwJEmSmm/WrHR1JliSJEm58eab6dqvX7Y5CmQJliRJUvPNmgVt28Lf/V3WSQpiCZYkSVLzzZqVRiFaVUa9rIyUkiRJKm+zZlXMKARYgiVJklQMlmBJkiTlypIlsHSpJViSJEk5Ur89miVYkiRJufHGG+laIXsEgyVYkiRJzfXKKxAC9O+fdZKCWYIlSZLUPK+8AnvuCR06ZJ2kYJZgSZIkNc/MmbDvvlmnaBJLsCRJkrbfhg3w6quwzz5ZJ2kSS7AkSZK23+zZsGqVK8GSJEnKkZkz09USLEmSpNywBEuSJCl3pk+Hbt1g112zTtIklmBJkiRtv8mT4YADsk7RZJZgSZIkbZ/aWpg2zRIsSZKkHHn1VVizBoYMyTpJk1mCJUmStH0mT05XV4IlSZKUG5MnQ9u2FXdQBliCJUmStL0mTYKBA6GmJuskTWYJliRJUtPFCOPHw4EHZp1ku1iCJUmS1HRvvAGLF8PBB2edZLtYgiVJktR0L76YrsOGZZtjO1mCJUmS1HQvvggdOqSZ4ApkCZYkSVLTvfgifPzj0KZN1km2iyVYkiRJTbNmDUycWLGjEGAJliRJUlONH5+K8GGHZZ1ku1mCJUmS1DRPPZWuhx+ebY5msARLkiSpaZ56CgYNgl13zTrJdrMES5IkqXC1tfDMM3DEEVknaRZLsCRJkgo3cSKsWAGf/GTWSZrFEixJkqTCjRuXrpZgSZIk5cZTT8GAAdCjR9ZJmsUSLEmSpMKsXZtK8Kc+lXWSZrMES5IkqTDPPgsffgjHHZd1kmazBEuSJKkwDz0ENTVw9NFZJ2k2S7AkSZIKM2ZMOiVup52yTtJslmBJkiRt29y5MGUKHH981kmKwhIsSZKkbRszJl2rYB4YLMGSJEkqxP/+L/TsmY5LrgKWYEmSJG3dihVpJfikkyCErNMUhSVYkiRJWzdmDKxaBaeemnWSoim4BIcQbgohLAghTGv03HdDCHNDCJPqHiNbJqYkSZIyc9dd0LUrHH541kmKpikrwTcDm5uEvjrGOKTu8WBxYkmSJKksrFkDDzyQRiHatMk6TdEUXIJjjOOAxS2YRZIkSeXmscfSKXGnnJJ1kqIqxkzwhSGEKXXjEp2L8PMkSZJULm69FXbZpSpOiWusuSX4WmAvYAgwH7hqS28MIYwOIYwPIYxfuHBhM3+tJEmSWtzy5XDPPXDGGdC2bdZpiqpZJTjG+F6McX2McQNwAzBsK++9PsY4NMY4tFu3bs35tZIkSSqFu++GlSth1KiskxRds0pwCGH3Rt+eDEzb0nslSZJUYW65Bfr2hU98IuskRVfwR/xCCLcCRwJdQwjvAP8OHBlCGAJE4C3gqy2QUZIkSaU2dy48/jhccUXVHJDRWMElOMb4hc08fWMRs0iSJKlc3HwzxFiVoxDgiXGSJEna1Pr1cP31aUeI/v2zTtMiLMGSJEna2EMPwezZcP75WSdpMZZgSZIkbezaa6FHDzjxxKyTtBhLsCRJkhq8+WZaCT73XKipyTpNi7EES5IkqcENN6TdIEaPzjpJi7IES5IkKVmxIn0g7oQTYI89sk7ToizBkiRJSn73O3j/ffjGN7JO0uIswZIkSYLaWrjqqnQ63KGHZp2mxRV8WIYkSZKq2J13wltvwc9+lnWSknAlWJIkKe9ihB//GPbZJ80D54ArwZIkSXk3Zgy8/DL89rfQKh9rpPn4TylJkqTNixG+/W3o0wdGjco6Tcm4EixJkpRn994LEyaknSHats06Tcm4EixJkpRXGzbAd74DAwbAWWdlnaakXAmWJEnKqzvugKlT4Y9/hDb5qoWuBEuSJOXRmjXwrW/BwIFwxhlZpym5fFV+SZIkJb/8JcyalXaGaN066zQl50qwJElS3ixcCD/4ARx/PBx7bNZpMmEJliRJypvvfhdWrICf/jTrJJmxBEuSJOXJ5Mlw3XXw1a/CfvtlnSYzlmBJkqS8WL8+ld8uXdI4RI75wThJkqS8uO46eOEFuOWWVIRzzJVgSZKkPJg/H775TTj6aPjiF7NOkzlLsCRJUh5cemnaG/jaayGErNNkzhIsSZJU7e65B26/Ha64Avr3zzpNWbAES5IkVbMFC9KH4Q48EP7lX7JOUzYswZIkSdUqRjjvPPjgg/RhuJqarBOVDXeHkCRJqlY33wz33Qf//d+53hN4c1wJliRJqkZvvAGXXAJHHpmu2oglWJIkqdqsXg2f/zy0aZNWg1tZ+TblOIQkSVK1+frXYeLENAqx555ZpylL/muBJElSNbn9drjmGrj8cjjhhKzTlC1LsCRJUrV49VU491wYPhx++MOs05Q1S7AkSVI1WLIkrfy2bw+33eZ2aNvgTLAkSVKlq62FM86At96CJ56A3r2zTlT2LMGSJEmV7utfh0cfhd/+Fg47LOs0FcFxCEmSpEp2443w85/DpZfCOedknaZiWIIlSZIq1ZgxcP758OlPw09+knWaimIJliRJqkQvvQSnnQaDBsEdd6SDMVQwS7AkSVKlee01GDkSuneHhx6CTp2yTlRxLMGSJEmV5N134dhj09cPPww9emSbp0K5bi5JklQpFi6EESNgwQJ48kno3z/rRBXLEixJklQJ3n8/FeA33oAHH4SDDso6UUWzBEuSJJW7JUvgmGPSscj33w+f+lTWiSqeJViSJKmcLV0Kxx0H06fDX/6SyrCazRIsSZJUrhYsSHsAz5gBd94Jxx+fdaKqYQmWJEkqR3PmpFXf2bPTCET9jhAqCkuwJElSuXn9dTj66DQK8cgjcNhhWSeqOpZgSZKkcjJlSlr1XbcubYN24IFZJ6pKBR+WEUK4KYSwIIQwrdFzXUIIj4YQXqu7dm6ZmJIkSTlQv+rbujWMG2cBbkFNOTHuZuC4TZ77V+DxGGN/4PG67yVJktRUN90En/kM9O0Lzz8P++2XdaKqVnAJjjGOAxZv8vSJwO/rvv49cFKRckmSJOVDjPDtb8M558BRR8HTT0OvXlmnqnpNWQnenN1ijPMB6q7dt/TGEMLoEML4EML4hQsXNvPXSpIkVYGVK+Gss+A//iOV4AcegE6dsk6VC80twQWLMV4fYxwaYxzarVu3Uv1aSZKk8vT222n+99Zb4Yc/hBtugJqarFPlRnN3h3gvhLB7jHF+CGF3YEExQkmSJFW1sWPh9NPTDhAPPAAjR2adKHeauxJ8H3B23ddnA/c28+dJkiRVrxjhl7+EESOga1d48UULcEaaskXarcBzwN4hhHdCCOcAVwLHhBBeA46p+16SJEmbWr4cvvQluPjitAvECy/AgAFZp8qtgschYoxf2MJLRxcpiyRJUnWaPBk+//l0Etz3vgdXXAGtSvbRLG2GJ8ZJkiS1lBjhuuvg0kuhSxd4/HE48sisU4kS7g4hSZKUK8uWwZlnwte+lorvpEkW4DJiCZYkSSq2J56A/feHu+6CK6+EBx+E7ls8TkEZsARLkiQVy6pVcNllcPTRsMMO8Mwz8C//4vxvGXImWJIkqRgmTIBRo2DmTLjwQvjRj6BDh6xTaQv81xJJkqTmWLs27fgwfDh88AE8/HDaC9gCXNZcCZYkSdpeL7wA554L06bBF74A11wDnTtnnUoFcCVYkiSpqZYvh0sugUMOgaVL4b774E9/sgBXEFeCJUmSmuKhh+D882HOHLjgAviv/4JOnbJOpSZyJViSJKkQc+bAGWfAyJGw447w9NPwq19ZgCuUJViSJGlr1qyBH/4Q9tknjT1873vw8stw6KFZJ1MzOA4hSZK0JWPGwMUXw2uvwUknwdVXQ58+WadSEbgSLEmStKk334STT4bjj0/fP/QQ3HOPBbiKWIIlSZLqLV4Ml1+eRh8eeSSNQUydCscdl3UyFZnjEJIkSWvWwK9/DT/4Qdry7Ctfge9/H3r2zDqZWogrwZIkKb9ihNtvh/32g3/+Zxg2DCZNghtvtABXOUuwJEnKpyefhE98Im171rFjOu54zBgYPDjrZCoBS7AkScqXZ5+Fo45Kjzlz4KabYOJE+PSns06mErIES5KkfJgwIR10ceihMGMG/Pzn8Prraf63deus06nELMGSJKm6TZ0Kp5wCQ4fCCy/Aj34Eb7yR9v9t3z7rdMqIu0NIkqTq9NJL8J//Cffem442/t734NJLPeZYgCVYkiRVm3HjUvl95BHYZRf4znfgkkugS5esk6mMWIIlSVLlizHt7vCf/wl//St07w5XXglf+5orv9osS7AkSapctbVw113wk5+kD7716gW/+AWccw506JB1OpUxS7AkSao8H3wAv/1t2uFh9mz4+7+HG26AL30J2rbNOp0qgCVYkiRVjjlzUvG94YZUhA8/PK38nnACtHLTKxXOEixJksrfhAlw1VXpiGOA006Dyy+Hgw7KNpcqliVYkiSVpzVr4M474Zpr4LnnYKed0i4PF18Me+6ZdTpVOEuwJEkqL3PmwHXXpZGHBQugf3+4+up0stvOO2edTlXCEixJkrIXIzz5ZFr1vfde2LABPvtZuPBCGDHCeV8VnSVYkiRlZ/Fi+J//gd/8BmbOhF13ha9/Hc4/H/r0yTqdqpglWJIkldaGDfDUU2mLs7vuSrO/Bx0EN98MZ5wB7dtnnVA5YAmWJEmlMX9+Kro33ghvvJGOND7vvHSwxZAhWadTzliCJUlSy6mthTFj0qrvAw/A+vVwxBHw3e/CqafCDjtknVA5ZQmWJEnFN2UK3HIL/PGPaQV4t93SrO8//iMMGJB1OskSLEmSiuTdd+FPf4I//AEmT4Y2beAzn4Gzz047PdTUZJ1Q+j+WYEmStP1Wrkxbmv3hD/DII+lDb8OGwa9+lT7k1rVr1gmlzbIES5Kkplm/HsaOTaMOd94JH34IvXvDN78JZ50F++yTdUJpmyzBkiRp2zZsgGefhT//ORXf995LxxiffjqMGgWf/KQHWqiiWIIlSdLmxQjjx8Ntt6XHO++kPXw/+1k480wYOdLdHVSxLMGSJKlBjDB1aiq9f/4zzJqVPtB23HFw5ZXwuc+lFWCpwlmCJUnKuxjh5Zfh7rvTCW6vvAKtW8NRR8G//RucfDJ07px1SqmoLMGSJOVR/Yzv3Xenx9tvp5neI46Aiy6C006D7t2zTim1GEuwJEl5sW5d2tXh7rvhL39J+/q2bQvHHAPf+U4adXBLM+WEJViSpGq2ahU8+mgac7j/fliyBHbcMX2o7ZRT0rVTp6xTSiVnCZYkqdrMmwcPPJBK7+OPpyLcuXNa6T3llLTy664OyjlLsCRJla7+g233358eEyak5/v0gXPPTeX3iCM8tlhqxBIsSVIlWrUqrfI+8EB6zJ0LIcAhh8APfwgnnAD77Zeek/QRRSnBIYS3gA+B9UBtjHFoMX6uJElq5O234eGHU+l97LFUhDt2hGOPTaV35Ejo1i3rlFJFKOZK8KdijIuK+PMkScq31ath3DgYMyY9Zs5Mz++5J5xzTiq+RxwB7dplm1OqQI5DSJJULmKE115rKL1jx6bV3nbtUtkdPTqd3Lb33o45SM1UrBIcgUdCCBG4LsZ4/aZvCCGMBkYD9O7du0i/VpKkCrd8OTzxREPxffPN9PyAAXDeean0HnEEdOiQbU6pyhSrBB8aY5wXQugOPBpCeCXGOK7xG+qK8fUAQ4cOjUX6vZIkVZb169PuDY89lh5//Ws6xGLHHeHoo+Eb30gzvv36ZZ1UqmpFKcExxnl11wUhhHuAYcC4rf9VkiTlQIzwt781lN6xY2Hp0vTa4MFw6aVptffQQ53tlUqo2SU4hLAj0CrG+GHd158Gvt/sZJIkVar589P2ZfXFd+7c9Pyee8Kpp8KIEXDUUdC9e7Y5pRwrxkrwbsA9IQ3otwH+FGMcU4SfK0lSZfjgA3jqqYbSO2NGer5LlzTiMGJEuvbr5wfapDLR7BIcY5wFHFCELJIkVYYPPkizvGPHpseECbBhQzqK+PDD4ctfTsX3gAOgVauMw0raHLdIkyRpW5Ytg6efTqu9Y8fCxImp9NbUwPDh8K1vpZXeQw5xrleqEJZgSZI2tXTpxqX35ZdT6W3bNpXeK65I25YNH+7WZVKFsgRLkvT++/DMM6nwPvVUKr0xplXd4cPh29+GI4+Egw9OIw+SKp4lWJKULzHCW2+lmd76R/0H2dq1SyMN//7vDaW3ffss00pqIZZgSVJ1W78epkzZuPTOm5de23ln+MQn4ItfhMMOg2HDLL1STliCJUnVZcUKePHFhsL73HPw4YfptT32SLO8hx2WHgMHQuvW2eaVlAlLsCSpss2bl4rus8+m0jtxItTWpv14Bw2Cs85qKL29e2edVlKZsARLkirHmjXpQ2vPP5+K7/PPw+zZ6bV27dI4wze+kQrvIYdA587Z5pVUtizBkqTyNWfOxoV3wgRYuza9tsceqehedlnaweFjH3OPXkkFswRLksrD6tVplKG+8D73HMydm15r1w6GDoWLL06Fd/hw6Nkz27ySKpolWJJUejHCm2+mD7DVF96XX4Z169LrffrAJz+Zyu4hh6Tjh9u2zTSypOpiCZYktbx334WXXkqPF1+E8ePTARWQDp846KA01nDIIan49uiRbV5JVc8SLEkqrg8+SCW3cemdMye91qpV2pbspJNS8T3oINh/f6ipyTazpNyxBEuStt+aNTB5ciq69YX31VfTuANAv37pMIphw1LhPfBA2HHHbDNLEpZgSVKh1q+HV17ZuPBOmdIwx7vbbqno/sM/pOvQodC1a7aZJWkLLMGSpI+qrYUZM9KWZBMnpuukSbBqVXp9p51Syb3ssoZV3j32SAdUSFIFsARLUt6tXQvTpzeU3QkT0grv6tXp9R13THvwjh6dxhkOOgj23jvN90pShbIES1KerFkDU6duXHinTm04gGKnnVLRveCCdP34x6F/f2jdOtvcklRklmBJqlarVqUV3frCO3EiTJvWMMO7yy6p6F5ySUPh3WsvV3gl5YIlWJKqwdKlaZeGSZPSY+LENOKwfn16vUuXVHIvv7yh8Pbt6wyvpNyyBEtSJYkRZs9uKLv1j7feanhPjx4wZAiccEJD4e3d28IrSY1YgiWpXK1dCzNnfrTwLl2aXg8BBgyAgw+G889PxfeAAzxtTZIKYAmWpHKw6TjDpElpnKF+fneHHWDwYDjjjFR2hwxJJ6158IQkbRdLsCSVUozpCOH6ovvyyx8dZ+jePW1JduyxDYXXHRokqagswZLUUlasSLsxTJmStiGbMiU9lixJrzceZ/jqVxsKr+MMktTiLMGS1FwbNsCsWQ0lt77wvvFGWvkF6NgxjS+cfvrG4wwdO2abXZJyyhIsSU2xePHGq7pTp6bHypXp9RDS6MKQITBqVJrjHTwY+vRx/11JKiOWYEnanHXr4NVXG8pufeF9552G93TpknZjOPfchrI7cCB06JBdbklSQSzBkvItRnj33Y3L7pQpaWuy+p0Zampg333hyCNT0d1//3TdfXf33pWkCmUJlpQfH34IM2akD6s1Hml4//2G9/TsmQru8cc3lN2994a2bbPLLUkqOkuwpOqzenUaZZg2beNH423IOnSAQYPg5JMbVnf33x923TWz2JKk0rEES6pctbVpB4ZNy+5rr8H69ek9bdrAPvvA8OFpdnfgwFR++/Z1311JyjFLsKTyt2EDzJ790bI7c2Y6WhjSbO5ee6WCe9pp6TpoUNqpwVEGSdImLMGSykeM8N57Hy2706fD8uUN79tjj1Rwjzmmoezuu6+7MkiSCmYJlpSNJUtSud208Db+kFrXrmlO9ytfaSi7++0Hu+ySXW5JUlWwBEtqWStWNOzI0Pgxb17Dezp1SgX31FMbZnYHDYLu3bPLLUmqapZgScWxdu3md2R4882Go4Pbt08ruSNGNBTdQYOgVy/325UklZQlWFLTrF8Ps2Z9tOz+7W9ptwZIOzLsvTcMHQpf/nJD2e3Xzx0ZJEllwRIsafPqd2Son9utv86cmfbhhbR6269fw3679WV3wAB3ZJAklTVLsJR3McL8+RvvxDBtWprjbbwjQ69eqeAeddTGOzLsuGN22SVJ2k6WYClPFi7cuOjWX5cubXjPbrulD6f94z82fEjNHRkkSVXGEixVo6VLP1p0p01LJbhe586p4J55ZroOHJge3bpll1uSpBKxBEuVbPnyNLbQuOhOnw5z5za8p2PHVG4/97mNtx/r0cMdGSRJuWUJlirBqlXwyisfXd19662G99RvP3b00Q1ld+BA6N3bsitJ0iYswVI5WbcubTW26dzu66+n3RoAampgn31g+HA499yGwtu3r9uPSZJUIEuwlIVN99qtL7yvvtqw126rVtC/fzo2uPHcbv/+qQhLkqTtZgmWWtp778HUqRs/pk9PIw71+vX76Nzu3nunEQdJklR0lmCpWFasSOV208LbeGGSvL8AAA5pSURBVEeG3XZLK7vnn5+u9duPudeuJEklVZQSHEI4Dvg50Br4bYzxymL8XKksrV+fZnQ3LbtvvJEOngDo0CGt6J5wQiq79Y/u3bPNLkmSgCKU4BBCa+Aa4BjgHeClEMJ9McYZzf3ZUqZibBhlmDKloezOmNFwbHD93O6QITBqVEPZ7dcvvSZJkspSMVaChwGvxxhnAYQQ/gycCFiCVTlWr06jDJMmbVx4Fy1qeE+PHqngXnBBQ9ndbz/YYYfsckuSpO1SjBLcE5jT6Pt3gIOL8HOL6o477uCdd97JOobKwdq16ZCJ5cvTHO/y5bByZcMoQ+vW0LMnDBiQZnU7dkzXxjsyLFkC48alhyRJ2qpevXpx+umnZx1jI8UowZvbhT9+5E0hjAZGA/Tu3bsIv1bahhhTuW1cdpcvTyW4Xvv2qeR27ZquHTum5zxcQpKkqlaMEvwOsEej73sB8zZ9U4zxeuB6gKFDh36kJLe0cvu3DxXZsmUweXLDY9KkNN5QP7vbtm0aXRgyBA44oOHRpUu2uSVJUiaKUYJfAvqHEPoCc4EzgX8ows+VNu+992DCBJg4MT1efnnj44O7dk1l95/+KRXdIUPSCWseMCFJkuo0uwTHGGtDCBcCD5O2SLspxji92cmkGGHu3IayW1985zX6g4YBA2DYMBg9umGVd/fdHWeQJElbVZR9gmOMDwIPFuNnKadihLff3niFd+JEWLAgvd6qVVrNPeoo+PjH4cADU+nt1Cnb3JIkqSJ5YpxKL0aYMwdefBFeeqmh+C5Zkl5v3TodNPGZz6Sye+CBaYXXU9UkSVKRWILV8pYsSWX3xRcbHu+9l16rqUn77Z52WkPhHTw47dAgSZLUQizBKq7Vq9PuDI0L79/+1vD6PvvAscemOd5hw1Lhbdcuu7ySJCmXLMHafjGmXRmeeQaeey4V3smTYd269Pruu8PBB8OXv5wK78c/DrvskmViSZIkwBKsplizJm1H9uyzqfg++yy8+256rWPHVHQvv7xhlbdnz2zzSpIkbYElWFu2cGEquvWPl15KRRigXz8YMQI+8Yn0GDQofaBNkiSpAliC1WDuXHjqqfQYO7ZhlremBoYOhQsvTIX3kEPSqIMkSVKFsgTn2Zw5DYX3qafg9dfT8zvvDIcfDuecA4cemmZ53a1BkiRVEUtwnsybB4891lB6Z81Kz++yC3zyk3DBBXDEEWlPXkcbJElSFbMEV7OVK2HcOHjkkfSYXneadZcuqfRefHEqvfvvb+mVJEm5YgmuJhs2wKRJqfA++ij89a+wdm3ah/fww+Hss+GYY9LevK1aZZ1WkiQpM5bgSrd8eSq8998P//u/sGBBen7w4LTSe8wxqQDvsEO2OSVJksqIJbgSvf12Kr0PPABPPplWe3feGY4/HkaOTMW3R4+sU0qSJJUtS3ClmDkT7rgD7rwTpk5Nzw0YABddBJ/9bNrFoaYm24ySJEkVwhJcrmJMH2S7885UfmfMgBBS2b3qqlR8BwzIOqUkSVJFsgSXm1mz4JZb4M9/hldeScX38MPhl7+EU06Bv/u7rBNKkiRVPEtwOVi2LK32/v73aUeHENIWZhddlIqv872SJElFZQnOSozpQ23XXw/33gurV8M++8B//Rd88YvQu3fWCSVJkqqWJbjUlixJK76/+Q28+ip07gznngtf+hIMHZpWgSVJktSiLMGlMm0aXH013HorrFoFw4enMnz66e7hK0mSVGKW4JYUIzz9NPz4x+kgix12gFGj4Pzz4WMfyzqdJElSblmCW0KM8NBD8P3vwwsvQNeu8L3vwQUXpK8lSZKUKUtwsT3xBFxxBTz3HPTtC7/+NZx9NnTokHUySZIk1bEEF8uUKXDZZakE9+oF110HX/mKp7hJkiSVoVZZB6h4ixfDhRemGd9Jk+BnP4PXXoPRoy3AkiRJZcqV4O0VI/zud/D//l/a9uxrX0szwF26ZJ1MkiRJ22AJ3h7z5sF558GDD6YjjX/1Kxg8OOtUkiRJKpDjEE11220wcGA67e0Xv4CxYy3AkiRJFcYSXKi1a+Gii+DMM2HffWHy5PR9K/8rlCRJqjSOQxRiwQI46aS07dlll8GPfuSH3iRJkiqYJXhbXn8djjsuzQHfdht8/vNZJ5IkSVIzWYK3ZuLEVIA3bEj7/w4fnnUiSZIkFYEleEumTIERI6BTJ3j4Ydh776wTSZIkqUgswZvzt7/BMceko46ffDIdfyxJkqSq4dYGm1q6FE44IY1APP64BViSJKkKuRLc2Pr1aQu0WbPSDLAjEJIkSVXJEtzYT3+a5n9/85t0EpwkSZKqkuMQ9aZNg+98B049FUaPzjqNJEmSWpAlGKC2Fs4+G3beGa69FkLIOpEkSZJakOMQADfemPYEvv126NYt6zSSJElqYa4Ef/hhGoM4/HA47bSs00iSJKkEXAm+9lpYsADuv98xCEmSpJzI90rwmjXws5+lk+GGDcs6jSRJkkok3yvBt94K8+fDzTdnnUSSJEkllO+V4N/9DgYMSEckS5IkKTfyW4LffhvGjYNRo5wFliRJypn8luA//Sldv/jFbHNIkiSp5PJbgu+5B4YPh759s04iSZKkEmtWCQ4hfDeEMDeEMKnuMbJYwVrU0qUwYQJ8+tNZJ5EkSVIGirE7xNUxxp8W4eeUzrhxsGEDHHVU1kkkSZKUgXyOQzzxBLRvn8YhJEmSlDvFKMEXhhCmhBBuCiF0LsLPa3lPPAGHHgrt2mWdRJIkSRnYZgkOITwWQpi2mceJwLXAXsAQYD5w1VZ+zugQwvgQwviFCxcW7T9Ak61YAdOmpRIsSZKkXNrmTHCMcUQhPyiEcAPwwFZ+zvXA9QBDhw6NhQYsuunTIUYYMiSzCJIkScpWc3eH2L3RtycD05oXpwQmT07XwYOzzSFJkqTMNHd3iB+HEIYAEXgL+GqzE7W06dOhQwf3B5YkScqxZpXgGOOoYgUpmVmzoF8/aJXPjTEkSZKUxy3S3nwzlWBJkiTlVr5KcIypBDsKIUmSlGv5KsGLFqUt0izBkiRJuZavEjxvXrr26pVtDkmSJGUqXyV40aJ07do12xySJEnKVL5K8Pvvp6slWJIkKdfyVYLrV4J33TXbHJIkScpUvkpw/UqwJViSJCnX8lWCFy2CTp2gpibrJJIkScpQvkrw++87DyxJkqScleBFixyFkCRJEm2yDlBSffvCXntlnUKSJEkZy1cJvvbarBNIkiSpDORrHEKSJEnCEixJkqQcsgRLkiQpdyzBkiRJyh1LsCRJknLHEixJkqTcsQRLkiQpdyzBkiRJyh1LsCRJknLHEixJkqTcsQRLkiQpdyzBkiRJyh1LsCRJknLHEixJkqTcsQRLkiQpdyzBkiRJyh1LsCRJknLHEixJkqTcCTHG0v/SEBYCb5f8F0NXYFEGv1el5X3OB+9zPnifq5/3OB+yvM97xhi7bfpkJiU4KyGE8THGoVnnUMvyPueD9zkfvM/Vz3ucD+V4nx2HkCRJUu5YgiVJkpQ7eSvB12cdQCXhfc4H73M+eJ+rn/c4H8ruPudqJliSJEmC/K0ES5IkSdVZgkMIx4UQXg0hvB5C+NfNvN4uhHBb3esvhBD6lD6lmquA+/zPIYQZIYQpIYTHQwh7ZpFTzbOt+9zofaeFEGIIoaw+faxtK+QehxA+X/f38/QQwp9KnVHNV8A/s3uHEJ4MIbxc98/tkVnk1PYLIdwUQlgQQpi2hddDCOEXdf8bmBJCOLDUGRuruhIcQmgNXAMcD+wHfCGEsN8mbzsHWBJj/HvgauBHpU2p5irwPr8MDI0xDgbuBH5c2pRqrgLvMyGEnYCLgRdKm1DNVcg9DiH0B74JHBpjHAhcWvKgapYC/16+Arg9xvgx4Ezg16VNqSK4GThuK68fD/Sve4wGri1Bpi2quhIMDANejzHOijGuBf4MnLjJe04Efl/39Z3A0SGEUMKMar5t3ucY45MxxpV13z4P9CpxRjVfIX8/A/yA9C85q0sZTkVRyD0+D7gmxrgEIMa4oMQZ1XyF3OcIdKr7emdgXgnzqQhijOOAxVt5y4nAH2LyPLBLCGH30qT7qGoswT2BOY2+f6fuuc2+J8ZYCywDdi1JOhVLIfe5sXOAh1o0kVrCNu9zCOFjwB4xxgdKGUxFU8jfywOAASGEZ0IIz4cQtrbSpPJUyH3+LnBWCOEd4EHgotJEUwk19f+7W1SbrH5xC9rciu6mW2AU8h6Vt4LvYQjhLGAocESLJlJL2Op9DiG0Io00fblUgVR0hfy93Ib0x6dHkv5E5+kQwqAY49IWzqbiKeQ+fwG4OcZ4VQjhEOCWuvu8oeXjqUTKqn9V40rwO8Aejb7vxUf/SOX/3hNCaEP6Y5etLd+r/BRynwkhjAD+DfhcjHFNibKpeLZ1n3cCBgFjQwhvAcOB+/xwXEUp9J/Z98YY18UY3wReJZViVY5C7vM5wO0AMcbngPZA15KkU6kU9P/dpVKNJfgloH8IoW8IoS1puP6+Td5zH3B23denAU9EN0yuNNu8z3V/TH4dqQA7Q1iZtnqfY4zLYoxdY4x9Yox9SLPfn4sxjs8mrrZDIf/M/gvwKYAQQlfSeMSskqZUcxVyn2cDRwOEEPYlleCFJU2plnYf8KW6XSKGA8tijPOzClN14xAxxtoQwoXAw0Br4KYY4/QQwveB8THG+4AbSX/M8jppBfjM7BJrexR4n38CdATuqPvc4+wY4+cyC60mK/A+q4IVeI8fBj4dQpgBrAe+EWN8P7vUaqoC7/PlwA0hhMtIf0T+ZReoKksI4VbS2FLXutnufwdqAGKMvyHNeo8EXgdWAl/JJmniiXGSJEnKnWoch5AkSZK2yhIsSZKk3LEES5IkKXcswZIkScodS7AkSZJyxxIsSZKk3LEES5IkKXcswZIkScqd/w8MdMkURdSNiQAAAABJRU5ErkJggg==\n",
"text/plain": [
"