22_succ_aprx.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Foundations of Computational Economics #22\n",
    "\n",
    "by Fedor Iskhakov, ANU\n",
    "\n",
    "<img src=\"_static/img/dag3logo.png\" style=\"width:256px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## Successive approximations (fixed point iterations)\n",
    "\n",
    "<img src=\"_static/img/lecture.png\" style=\"width:64px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "<img src=\"_static/img/youtube.png\" style=\"width:65px;\">\n",
    "\n",
    "[https://youtu.be/AQt9Q9qc3io](https://youtu.be/AQt9Q9qc3io)\n",
    "\n",
    "Description: Scalar and multivariate solver. Equilibrium in market of platforms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### What is successive approximations?\n",
    "\n",
    "- root finding method for the equations of the type  \n",
    "\n",
    "\n",
    "$$\n",
    "x = F(x) \\quad \\Leftrightarrow \\quad F(x) - x = 0\n",
    "$$\n",
    "\n",
    "- **fixed point equation**  \n",
    "- conversion is always possible (although does not guaranteed to work as we’ll see below)  \n",
    "\n",
    "\n",
    "$$\n",
    "g(x) = 0 \\quad \\Leftrightarrow \\quad \\underbrace{g(x) + x}_{F(x)} = x\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Main algorithm\n",
    "\n",
    "$$\n",
    "x = F(x) \\rightarrow x_{n+1} = F(x_n)\n",
    "$$\n",
    "\n",
    "1. Initialize the iterations at $ x_0 $  \n",
    "1. Iterate on $ x_{i+1} = F(x_i) $ until $ x_{i+1} $ is close enough to $ x_i $  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Stationary distribution of Markov chains\n",
    "\n",
    "For aperiodic and irreducible Markov chain with transition probability matrix $ P $, the stationary distribution\n",
    "is given by $ \\psi^\\star = \\psi^\\star P $\n",
    "\n",
    "Successive approximation applicable with $ F(\\psi) = \\psi P $\n",
    "\n",
    "Have seen this is last two videos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Will it converge?\n",
    "\n",
    "*The central question for the successive approximations method*\n",
    "\n",
    "- it will from any starting point for aperiodic and irreducible Markov chains  \n",
    "- but what are the general conditions?  \n",
    "\n",
    "\n",
    "Will come back to this question in general terms when talking about the dynamic models and Bellman operator\n",
    "(which is a *contraction mapping*)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Scalar example\n",
    "\n",
    "$$\n",
    "\\frac{1}{2} - \\exp\\big(-(x-2)^2\\big) = 0\n",
    "$$\n",
    "\n",
    "Using the trick above convert to $ x = F(x) $\n",
    "\n",
    "$$\n",
    "F(x) = x - \\exp\\big(-(x-2)^2\\big) + \\frac{1}{2}\n",
    "$$\n",
    "\n",
    "$$\n",
    "x_{i+1} = x_i - \\exp\\big(-(x_i-2)^2\\big) + \\frac{1}{2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "eq = lambda x:  - np.exp(-(x-2)**2) + .5  # initial equation\n",
    "F = lambda x: x - np.exp(-(x-2)**2) + .5  # fixed point equation\n",
    "dF = lambda x: 1 + np.exp(-(x-2)**2)*2*(x-2)  # derivative (for later)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def solve_sa(F,x0,tol=1e-6,maxiter=100,callback=None):\n",
    "    '''Computes the solution of fixed point equation x = F(x)\n",
    "    with given initial value x0 and algorithm parameters\n",
    "    Method: successive approximations\n",
    "    '''\n",
    "    for i in range(maxiter):  # main loop\n",
    "        x1 = F(x0)  # update approximation\n",
    "        err = np.amax(np.abs(x0-x1))  # allow for x to be array\n",
    "        if callback != None: callback(iter=i,err=err,x=x1,x0=x0)\n",
    "        if err<tol:\n",
    "            break  # break out if converged\n",
    "        x0 = x1  # get ready to the next iteration\n",
    "    else:\n",
    "        raise RuntimeError('Failed to converge in %d iterations'%maxiter)\n",
    "    return x1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter 0, x = 1.1321205588 (err = 1.321e-01)\n",
      "iter 1, x = 1.1612700962 (err = 2.915e-02)\n",
      "iter 2, x = 1.1664038940 (err = 5.134e-03)\n",
      "iter 3, x = 1.1672707853 (err = 8.669e-04)\n",
      "iter 4, x = 1.1674161464 (err = 1.454e-04)\n",
      "iter 5, x = 1.1674404922 (err = 2.435e-05)\n",
      "iter 6, x = 1.1674445689 (err = 4.077e-06)\n",
      "iter 7, x = 1.1674452515 (err = 6.826e-07)\n",
      "iter 8, x = 1.1674453659 (err = 1.143e-07)\n",
      "iter 9, x = 1.1674453850 (err = 1.914e-08)\n",
      "iter 10, x = 1.1674453882 (err = 3.205e-09)\n",
      "iter 11, x = 1.1674453887 (err = 5.366e-10)\n",
      "iter 12, x = 1.1674453888 (err = 8.986e-11)\n",
      "Residual in the original equation = 1.50464e-11\n"
     ]
    }
   ],
   "source": [
    "def printme(**kwargs):\n",
    "    print('iter %d, x = %1.10f (err = %1.3e)'%(kwargs['iter'],kwargs['x'],kwargs['err']))\n",
    "x = solve_sa(F,x0=1.0,tol=1e-10,callback=printme)\n",
    "print('Residual in the original equation = %1.5e' % eq(x))\n",
    "# try different starting points, including >2.8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# plot to understand what is going on!\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.rcParams['figure.figsize'] = [12, 8]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAHSCAYAAAAezFYoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeVhV16H+8e8GQRRUnGec4zzjEOdZUSBJo3GIGUzatLlp7m16O7e/pE2H2za97U1u2zRpmqZt2qRtmuEcREBxFiecxzgPCE4oiMp81u+PhVdiNKIe3OfA+3me/RzO3ht4aam+rq69lmOMQURERERErgpxO4CIiIiISKBRSRYRERERuYZKsoiIiIjINVSSRURERESuoZIsIiIiInINlWQRERERkWvUcTvA9TRr1sx07NjR7RgiIiIiUoNt2rTprDGm+fWuBWRJ7tixI5mZmW7HEBEREZEazHGcoze6pukWIiIiIiLXUEkWEREREbmGSrKIiIiIyDVUkkVERERErqGSLCIiIiJyDZVkEREREZFrqCSLiIiIiFxDJVlERERE5BoqySIiIiIi11BJFhERERG5hkqyiIiIiMg1VJJFRERERK6hkiwiIiIicg2VZBERERGRa9y0JDuO095xnGWO4+xxHGeX4zj/cZ17HMdxXnEc54DjONsdxxlU6dpjjuPsrzge8/cPICIiIiLib3WqcE8Z8J/GmM2O4zQANjmOs9gYs7vSPXFAt4pjGPAqMMxxnCbAC0AsYCo+12OMOe/Xn0JERERExI9uOpJsjMkxxmyu+LgA2AO0vea2+4A/G2sdEO04TmtgKrDYGHOuohgvBqb59ScQERERkaB14cIFysrK3I7xKbc0J9lxnI7AQGD9NZfaAscrvc+qOHej8yIiIiJSixljyMzM5De/+Q2rVq1yO86nVGW6BQCO40QB/wK+Yoy5cO3l63yK+Yzz1/v6TwFPAcTExFQ1loiIiIgEmdzcXLxeL0ePHqVz584MGDDA7UifUqWS7DhOGLYg/9UY8/51bskC2ld63w7Irjg/7przy6/3PYwxrwOvA8TGxl63SIuIiIhI8PL5fKxdu5bly5cTGhpKYmIiAwYMwHGuN67qrpuWZMem/gOwxxjzyxvc5gG+7DjOu9gH9/KNMTmO46QCP3Ecp3HFfVOAb/sht4iIiIgEkZMnT+LxeMjJyaFHjx5Mnz6dBg0auB3rhqoykjwSeATY4TjO1opz3wFiAIwxvwOSgenAAeAysKDi2jnHcX4IbKz4vBeNMef8F19EREREAllZWRkrV65kzZo11KtXj5kzZ9KrV6+AHD2u7KYl2RizmuvPLa58jwGeucG1N4E3byudiIiIiASt48eP4/F4OHv2LP3792fKlCnUr1/f7VhVUuUH90REREREqqKkpISlS5eyfv16GjZsyLx58+jWrZvbsW6JSrKIiIiI+M2hQ4fwer3k5eUxZMgQJk6cSN26dd2OdctUkkVERETkjhUWFpKWlsbWrVtp2rQpjz/+OB06dHA71m1TSRYRERGRO7Jnzx6Sk5O5dOkSI0eOZNy4cdSpE9w1M7jTi4iIiIhrLl68yKJFi9i9ezetWrVi3rx5tG7d2u1YfqGSLCIiIiK3xBjD9u3bSUlJobS0lAkTJjBixAhCQ0PdjuY3KskiIiIiUmX5+fkkJSVx4MAB2rdvT2JiIs2aNXM7lt+pJIuIiIjITRlj2LhxI+np6RhjiIuLY8iQIQG/KcjtUkkWERERkc909uxZvF4vx44do0uXLsTHxxMdHe12rGqlkiwiIiIi11VeXs7atWtZvnw5YWFh3HffffTv37/Gjh5XppIsIiIiIp+Sk5ODx+Ph5MmT9OzZk+nTpxMVFeV2rLtGJVlERERE/k9ZWRkrVqxgzZo11K9fn1mzZtGrVy+3Y911KskiIiIiAsCxY8fweDzk5uYyYMAApkyZQr169dyO5QqVZBEREZFarqSkhPT0dDZs2ECjRo2YP38+Xbp0cTuWq1SSRURERGqxAwcOkJSURH5+PkOHDmXixImEh4e7Hct1KskiIiIitVBhYSGpqals27aNZs2asWDBAmJiYtyOFTBUkkVERERqmd27d5OcnMzly5cZPXo0Y8aMoU4d1cLK9J+GiIiISC1RUFDAokWL2LNnD61atWL+/Pm0atXK7VgBSSVZREREpIYzxrBt2zZSU1MpLS1l4sSJjBgxgpCQELejBSyVZBEREZEaLC8vD6/Xy6FDh4iJiSExMZGmTZu6HSvgqSSLiIiI1EA+n4+NGzeSnp6O4zhMnz6d2NjYWrGltD+oJIuIiIjUMGfOnMHr9XL8+HG6du1KfHw8jRo1cjtWUFFJFhEREakhysvLWbNmDStXriQ8PJz777+ffv36afT4Nqgki4iIiNQAOTk5fPTRR5w6dYrevXszbdo0oqKi3I4VtFSSRURERIJYaWkpK1asICMjg8jISGbPnk2PHj3cjhX0VJJFREREgtTRo0fxer3k5uYycOBApkyZQkREhNuxagSVZBEREZEgU1xczJIlS8jMzCQ6OppHHnmEzp07ux2rRlFJFhEREQki+/fvJykpiQsXLjBs2DAmTJhAeHi427FqHJVkERERkSBw+fJlUlNT2b59O82bN+fJJ5+kXbt2bseqsVSSRURERAKYMYbdu3eTnJxMUVERY8aMYfTo0dSpoxpXnfSfroiIiEiAKigoIDk5mb1799KmTRsSExNp2bKl27FqBZVkERERkQBjjGHLli2kpaVRXl7O5MmTGT58OCEhIW5HqzVUkkVEREQCyPnz5/F6vRw+fJgOHTqQmJhIkyZN3I5V66gki4iIiAQAn8/Hhg0bWLp0KY7jMGPGDAYPHqwtpV2ikiwiIiListOnT+PxeDhx4gTdunUjPj6ehg0buh2rVlNJFhEREXFJeXk5q1evZuXKldStW5fPfe5z9OnTR6PHAUAlWURERMQFJ06cwOPxcPr0afr06cO0adOIjIx0O5ZUUEkWERERuYtKS0tZtmwZ69atIyoqijlz5tC9e3e3Y8k1VJJFRERE7pIjR47g9Xo5d+4cgwYNYvLkyURERLgdS65DJVlERESkmhUVFbFkyRI2bdpE48aNefTRR+nUqZPbseQz3LQkO47zJhAPnDbG9LnO9a8DD1f6ej2B5saYc47jHAEKgHKgzBgT66/gIiIiIsFg3759JCUlcfHiRe69917Gjx9PWFiY27HkJqoykvwW8Gvgz9e7aIx5CXgJwHGcBOA5Y8y5SreMN8acvcOcIiIiIkHl0qVLpKamsmPHDlq0aMHs2bNp27at27Gkim5ako0xKx3H6VjFrzcXeOdOAomIiIgEM2MMO3fuJCUlhaKiIsaNG8eoUaMIDQ11O5rcAr/NSXYcpz4wDfhypdMGSHMcxwCvGWNe99f3ExEREQk0Fy5cYOHChezbt4+2bduSmJhIixYt3I4lt8GfD+4lAGuumWox0hiT7ThOC2Cx4zh7jTErr/fJjuM8BTwFEBMT48dYIiIiItXLGMPmzZtZvHgx5eXlTJkyhWHDhhESEuJ2NLlN/izJc7hmqoUxJrvi9bTjOB8AQ4HrluSKUebXAWJjY40fc4mIiIhUm3PnzuH1ejly5AgdO3YkISGBJk2auB1L7pBfSrLjOI2AscD8SucigRBjTEHFx1OAF/3x/URERETc5vP5WLduHcuWLSM0NJSEhAQGDhyoLaVriKosAfcOMA5o5jhOFvACEAZgjPldxW0PAGnGmEuVPrUl8EHFL0od4G/GmBT/RRcRERFxx+nTp/noo4/Izs6me/fuTJ8+nYYNG7odS/yoKqtbzK3CPW9hl4qrfO4Q0P92g4mIiIgEmrKyMlavXs2qVauIiIjgwQcfpHfv3ho9roG0456IiIhIFWRlZeHxeDhz5gz9+vVj6tSp1K9f3+1YUk1UkkVEREQ+Q0lJCcuWLWPdunU0bNiQuXPncs8997gdS6qZSrKIiIjIDRw+fBiv18v58+eJjY1l0qRJ1K1b1+1YcheoJIuIiIhco6ioiLS0NLZs2UKTJk147LHH6Nixo9ux5C5SSRYRERGp5OOPP2bhwoVcvHiRESNGMG7cOMLCwtyOJXeZSrKIiIgIcOnSJRYtWsSuXbto2bIlc+bMoU2bNm7HEpeoJIuIiEitZoxhx44dpKSkUFJSwvjx4xk5ciShoaFuRxMXqSSLiIhIrZWfn8/ChQvZv38/7dq1IzExkebNm7sdSwKASrKIiIjUOsYYMjMzWbJkCcYYpk6dytChQwkJCXE7mgQIlWQRERGpVXJzc/F6vRw9epTOnTsTHx9P48aN3Y4lAUYlWURERGoFn8/H2rVrWb58OXXq1CExMZEBAwZoS2m5LpVkERERqfFOnjyJx+MhJyeHHj16MH36dBo0aOB2LAlgKskiIiJSY5WVlbFy5UrWrFlDvXr1mDVrFj179tTosdyUSrKIiIjUSMePH8fj8XD27Fn69+/PlClTqF+/vtuxJEioJIuIiEiNUlJSwtKlS1m/fj2NGjXi4YcfpmvXrm7HkiCjkiwiIiI1xsGDB0lKSiIvL48hQ4YwceJE6tat63YsCUIqySIiIhL0CgsLSUtLY+vWrTRt2pQFCxYQExPjdiwJYirJIiIiEtT27NlDcnIyly5dYtSoUYwdO5Y6dVRx5M7oN0hERESC0sWLF1m0aBG7d++mVatWzJs3j9atW7sdS2oIlWQREREJKsYYtm/fTkpKCqWlpUyYMIERI0YQGhrqdjSpQVSSRUREJGjk5eWRlJTEwYMHad++PYmJiTRr1sztWFIDqSSLiIhIwDPGsHHjRtLT0zHGEBcXx5AhQ7QpiFQblWQREREJaGfPnsXr9XLs2DG6dOlCfHw80dHRbseSGk4lWURERAJSeXk5GRkZrFixgrCwMO677z769++v0WO5K1SSRUREJODk5OTg8Xg4efIkvXr1Ii4ujqioKLdjSS2ikiwiIiIBo6ysjBUrVrBmzRrq16/PQw89RM+ePd2OJbWQSrKIiIgEhGPHjuHxeMjNzWXAgAFMmTKFevXquR1LaimVZBEREXFVcXEx6enpbNy4kejoaObPn0+XLl3cjiW1nEqyiIiIuObAgQMkJSWRn5/P0KFDmThxIuHh4W7HElFJFhERkbuvsLCQ1NRUtm3bRrNmzXjiiSdo376927FE/o9KsoiIiNxVu3fvJjk5mcLCQkaPHs2YMWOoU0eVRAKLfiNFRETkrigoKGDRokXs2bOH1q1bM3/+fFq1auV2LJHrUkkWERGRamWMYevWraSlpVFaWsrEiRMZMWIEISEhbkcTuSGVZBEREak2eXl5eL1eDh06RExMDImJiTRt2tTtWCI3pZIsIiIifufz+di4cSPp6ek4jsP06dOJjY3VltISNFSSRURExK/OnDmDx+MhKyuLrl27Eh8fT6NGjdyOJXJLVJJFRETEL8rLy1mzZg0rV64kPDycBx54gL59+2r0WIKSSrKIiIjcsezsbDweD6dOnaJ3797ExcURGRnpdiyR26aSLCIiIrettLSUFStWkJGRQWRkJLNnz6ZHjx5uxxK5YyrJIiIicluOHj2Kx+Ph3LlzDBw4kClTphAREeF2LBG/UEkWERGRW1JcXMySJUvIzMwkOjqaRx55hM6dO7sdS8SvblqSHcd5E4gHThtj+lzn+jjgI+Bwxan3jTEvVlybBrwMhAJvGGN+6qfcIiIi4oL9+/eTlJTEhQsXGD58OOPHjyc8PNztWCJ+V5WR5LeAXwN//ox7Vhlj4iufcBwnFPgNMBnIAjY6juMxxuy+zawiIiLiksuXL5Oamsr27dtp3rw5Tz75JO3atXM7lki1uWlJNsasdByn42187aHAAWPMIQDHcd4F7gNUkkVERIKEMYbdu3eTnJxMUVERY8aMYfTo0dSpoxmb4idr10L37tCkidtJPsFfm6bf6zjONsdxFjmO07viXFvgeKV7sirOiYiISBAoKCjg73//O++99x7R0dE89dRTjB8/XgVZ/OPCBXjmGRg5En7yE7fTfIo/fss3Ax2MMRcdx5kOfAh0A663cri50RdxHOcp4CmAmJgYP8QSERGR22GMYcuWLaSlpVFeXs7kyZMZPnw4ISH+GluTWu+jj2xBzs6GZ5+FF15wO9Gn3HFJNsZcqPRxsuM4v3Ucpxl25Lh9pVvbAdmf8XVeB14HiI2NvWGZFhERkepz/vx5vF4vhw8fpkOHDiQmJtIkwP5vcAliOTm2FP/rX9C3L7z/Pgwd6naq67rjkuw4TivglDHGOI4zFDuFIxfIA7o5jtMJOAHMAebd6fcTERER//P5fKxfv56lS5cSEhJCfHw8gwYN0pbS4h8+H7zxBnzjG1BUZKdXfO1rEBbmdrIbqsoScO8A44BmjuNkAS8AYQDGmN8BM4GnHccpAwqBOcYYA5Q5jvNlIBW7BNybxphd1fJTiIiIyG07ffo0Ho+HEydO0K1bN+Lj42nYsKHbsaSm2LoVnn4a1q2D8ePhtdegWze3U92UY/tsYImNjTWZmZluxxAREanRysvLWb16NStXriQiIoJp06bRp08fjR6Lf1y4AM8/D//7v9C0Kfz85/DYYxBAv1+O42wyxsRe75oeTxUREamFTpw4gcfj4fTp0/Tt25epU6cSGRnpdiypCYyBf/wDnnsOTp6EL34RfvzjgFvi7WZUkkVERGqR0tJSli1bxrp164iKimLOnDl0797d7VhSU+zdC//+77B4MQwaBB9+GLAP5t2MSrKIiEgtceTIETweD+fPn2fw4MFMmjSJiIgIt2NJTZCfDy++CK+8ApGR8Otfw5e+BKGhbie7bSrJIiIiNVxRURGLFy9m8+bNNG7cmEcffZROnTq5HUtqAp8P3noLvv1tOHMGnnzSTq1o0cLtZHdMJVlERKQG27dvH0lJSVy8eJF7772X8ePHExbAy25JEFm3zk6t2LgRRoyA5GQYPNjtVH6jkiwiIlIDXbp0iZSUFHbu3EmLFi2YPXs2bdu2dTuW1ATHjsF3vgN//Su0aQNvvw3z5gXUqhX+oJIsIiJSgxhj2LlzJykpKRQVFTFu3DhGjRpFaBDPDZUAceEC/PSn8Ktf2fff+Y6dZhEV5W6uaqKSLCIiUkNcuHCBhQsXsm/fPtq2bUtiYiItasDcUHFZWRn8/vfwwgt23vH8+XbecUyM28mqlUqyiIhIkDPGsGnTJhYvXozP52PKlCkMGzaMkJAQt6NJMDMGkpLgm9+EPXtgzBg77zj2untv1DgqySIiIkHs3LlzeL1ejhw5QqdOnUhISKBx48Zux5Jgt3o1fOtbsGaN3UL6ww8hMbHGzTv+LCrJIiIiQcjn87Fu3TqWLVtGaGgoCQkJDBw4UFtKy53Zvt3ONV64EFq3htdegwULoBauiKKSLCIiEmROnTqFx+MhOzub7t27M336dBo2bOh2LAlmhw/D88/bFSsaNbIP6D37LNSv73Yy16gki4iIBImysjJWrVrF6tWriYiIYObMmfTq1Uujx3L7jh61D+H98Y92tPib34RvfAM0ZUclWUREJBhkZWXh8Xg4c+YM/fr1Y+rUqdSvxaN8coeysuAnP4E33rDzjL/0JbucW5s2bicLGCrJIiIiAaykpIRly5axbt06GjZsyLx58+jWrZvbsSRYZWfbqRSvvWZXr/j85205bt/e7WQBRyVZREQkQB06dAiv10teXh6xsbFMmjSJunXruh1LgtGRI/Dzn8Obb0J5uX0Y77vfhQ4d3E4WsFSSRUREAkxRURFpaWls2bKFJk2a8Pjjj9NBZUZux8cfw3/9l30gz3Hg8cft0m6dO7udLOCpJIuIiASQvXv3snDhQi5dusTIkSMZO3YsYbVw+S25Q9u22TnH//wnRETAM8/A174G7dq5nSxoqCSLiIgEgEuXLrFo0SJ27dpFy5YtmTt3Lm30EJXcCmNg6VJ46SVITYUGDexqFc89B9qe/JapJIuIiLjIGMOOHTtISUmhpKSE8ePHM3LkSEJDQ92OJsGirMyOGL/0EmzZAi1b2mXdnn5aS7ndAZVkERERl+Tn55OUlMSBAwdo164diYmJNG/e3O1YEiwKCuyDeL/6lV3vuHt3+P3vYf58O8VC7ohKsoiIyF1mjCEzM5MlS5ZgjGHatGkMGTKEkJAQt6NJMDh4EP73f21BLiiAUaPglVcgPh70O+Q3KskiIiJ3UW5uLh6Ph2PHjtG5c2fi4+NprP9LXG7mynzjl1+GpCQIDYWHHoL/+A8YOtTtdDWSSrKIiMhd4PP5WLt2LcuXL6dOnTokJiYyYMAAbSktn62gAN5+G377W9i5E5o3h+99z+6Qpwc7q5VKsoiISDU7efIkHo+HnJwcevTowfTp02nQoIHbsSSQbdsGr75q1ze+eBEGDoQ//hHmzNF847tEJVlERKSalJWVsXLlStasWUO9evWYNWsWvXr1cjuWBKqiIrtKxauvwtq1tgzPmWNHjYcOtZuByF2jkiwiIlINjh8/jsfj4ezZs/Tv35+pU6dSr149t2NJINq8Gf7wB/jb3yAvD+65B375S3jsMWjSxO10tZZKsoiIiB+VlJSQnp7Ohg0baNSoEQ8//DBdu3Z1O5YEmvPn7VSKP/wBtm6FunXhwQfhySdh/HiNGgcAlWQRERE/OXjwIElJSeTl5TFkyBAmTpxI3bp13Y4lgaK0FNLS4C9/gQ8/hOJiO9f417+GefO08UeAUUkWERG5Q4WFhaSlpbF161aaNm3KggULiImJcTuWBAJjYNMmW4zfeQfOnIGmTeHzn7ejxgMHup1QbkAlWURE5A7s2bOH5ORkLl26xKhRoxg7dix16uiv11rvwAH4+9/t8m1790J4OCQkwCOPQFycfS8BTf8rFhERuQ0XL14kOTmZPXv20KpVK+bNm0fr1q3djiVuOnYM/vEPePddO3oMdje8116DWbM0nSLIqCSLiIjcAmMM27ZtIzU1ldLSUiZMmMCIESMIDQ11O5q44fhxeP99W44zMuy52Fh46SW7I56m3QQtlWQREZEqysvLIykpiYMHD9K+fXsSExNp1qyZ27Hkbtu3zxbj99+HjRvtub594cc/htmzoUsXd/OJX6gki4iI3IQxho0bN7JkyRIA4uLiGDJkiLaUri18PruWsccDH3xgt4cGu8HHf/0XPPAAdO/ubkbxO5VkERGRz3D27Fk8Hg/Hjx+nS5cuxMfHEx0d7XYsqW6XL8OSJeD1wsKFkJMDISEwejS88grcfz+0b+92SqlGKskiIiLXUV5eTkZGBitWrCAsLIz777+ffv36afS4JjtwAFJTITkZli6120Q3aADTptmVKeLiQNNrag2VZBERkWvk5OTg8Xg4efIkvXr1Ii4ujqioKLdjib9dvAjLl0NKij0OHrTnO3eGp56yxXjMGC3XVkupJIuIiFQoKytj+fLlZGRkEBkZyUMPPUTPnj3djiX+UlZmH7RLT7dTKdauhZISqF/fbgX9la/YUWNtIy6oJIuIiABw7NgxPB4Pubm5DBgwgClTplCvXj23Y8md8Plgxw5YtswW4xUroKDAXhswAP7932HqVLuWcUSEu1kl4Ny0JDuO8yYQD5w2xvS5zvWHgW9WvL0IPG2M2VZx7QhQAJQDZcaYWD/lFhER8Yvi4mLS09PZuHEj0dHRzJ8/ny5awis4lZXBli2wcqUtxKtWQV6evda1K8ybBxMn2lFjzS2Wm6jKSPJbwK+BP9/g+mFgrDHmvOM4ccDrwLBK18cbY87eUUoREZFqcODAAZKSksjPz2fYsGFMmDCBcM0/DR4XLsD69XYTj4wMO33iykhxt24wc6adUzx2rDb1kFt205JsjFnpOE7Hz7ieUentOqDdnccSERGpPpcvXyYtLY1t27bRrFkznnjiCdprOa/A5vPB/v2wYYMtwxkZdiqFzweOYzfzePhhW4jHjIE2bdxOLEHO33OSnwQWVXpvgDTHcQzwmjHmdT9/PxERkSozxrBnzx6Sk5MpLCxk9OjRjBkzhjp19IhOwMnJsYX4yrFxI+Tn22sNGsDw4fD88zBiBAwbBg0buptXahy//angOM54bEkeVen0SGNMtuM4LYDFjuPsNcasvMHnPwU8BRCj/0tERET8rKCggOTkZPbu3Uvr1q2ZP38+rVq1cjuWGANZWbBpk93VbvNm+/HJk/Z6aCj06wdz5tgd7oYOhZ497XmRauSXkuw4Tj/gDSDOGJN75bwxJrvi9bTjOB8AQ4HrluSKUebXAWJjY40/comIiBhj2Lp1K2lpaZSVlTFp0iTuvfdeQkJC3I5W+xQXw+7dsH07bNt29fVsxaNLISHQqxdMmQKDBsGQITBwIGiVEXHBHZdkx3FigPeBR4wx+yqdjwRCjDEFFR9PAV680+8nIiJSVefPnycpKYlDhw7RoUMHEhISaNq0qduxar7ycrsxx65d9ti5077u3WtXoABbfPv0gfvus4V40CA7Yly/vrvZRSpUZQm4d4BxQDPHcbKAF4AwAGPM74DngabAbyu26ryy1FtL4IOKc3WAvxljUqrhZxAREfkEn8/Hxo0bSU9Px3EcZsyYweDBg7WltL8VFdmH6fbutceePfbYu9deu6JTJ+jdGxIToX9/e3TtqikTEtAcYwJvZkNsbKzJzMx0O4aIiAShM2fO4PF4yMrKomvXrsTHx9OoUSO3YwWv8nI4fhz27bOF+Mrrxx/D4cN2TvEVHTtCjx62EPfpY1979gRt6S0BynGcTTfax0OP84qISI1QXl7OmjVrWLlyJeHh4TzwwAP07dtXo8dVUVwMR47YKRLXO0pKrt4bGWnXII6NhfnzbQnu0QPuuUdTJaRGUUkWEZGgl52djcfj4dSpU/Tu3Zu4uDgiIyPdjhU4iovtChJHj9oyfPjwJ1+zsz85IhwZCZ07Q/fuMGOGLcDdutnX1q3tusQiNZxKsoiIBK3S0lKWL1/O2rVriYyMZPbs2fTo0cPtWHeXzwdnztgpEceP2zJ8/LgtxMeO2deTJz9ZgkNCoH17Oz1i8mT72rkzdOlijxYtVISl1lNJFhGRoHT06FE8Hg/nzp1j0KBBTJ48mYiICLdj+VdJiYMnK8IAACAASURBVC24J07YIyvr0x9nZX1yOgRAeLjdhrlDB5g2zb526GDPdeoE7dpBWJg7P5NIkFBJFhGRoFJcXMySJUvIzMykcePGPProo3Tq1MntWLfmSvnNyfn0ceKEnf6QnW1HiK9Vty60bWuL7rBhMHOm/bh9+6tH8+YaCRa5QyrJIiISNPbv309SUhIFBQUMHz6c8ePHEx4e7nYsyxg4f96W31Onrpbgyq9XPs7N/fTnO46d5tC2rR3xHT4c2rS5elwpxk2aqACL3AUqySIiEvAuX75MSkoKO3bsoHnz5syaNYt27drdnW9+8eIni++V43rvr532AHbkt1Ur+8Bb164wevTV95WPFi2gjv5aFgkU+l+jiIgELGMMu3btYtGiRRQVFTF27FhGjRpFnTstk8XFttReW3SvPU6dgkuXPv35ISF2SkOrVtCypV0GrVWrTx4tW9ryGx2tkV+RIKSSLCIiAamgoICFCxfy8ccf06ZNGxITE2nZsuWNP8EYuHDh+vN8K091OHkSzp27/tdo2tSW21at7HSHK2W38murVtCsmXaLE6nhVJJFRCSgGGPYsmULaWlplJeXM3nyZIb37UtITg7s3n31obYrR07O1Y8LCz/9BSMi7Ihuq1Z204tx4+z7KyO9V4pvixZ2VQgREVSSRUTETRUPuw388wjOXD5LJ180HXIM/5nRkEdLS2lRVkadl16C/PxPf25k5NWH2oYO/eT83jZtrn7csKGmO4jILVNJFhGR6lNUBMePM/DDOM4U5dLV1wiKiqG4qOK1GHw+dlc8g9c1K5emJ6H3ieaEdeqE07at3eyibdtPHm3aQIMG7v5sIlKjqSSLiMjtKyqyO7odPnx1i+MjR+y5Kzu9AWe+ChfDgZN5dkpD3boQFQVNm1Jepw4Ny09TrzyCrwz+ARPmzSO8SRM3fyoREZVkERH5DMbYDS0OHrTHoUNXXw8fthtfVFZ5p7cZM/5vp7eu534JEREs/8HK/5v3W15ezqpVq1i1ahWJ9SOYNm0affr0wdHUCBEJACrJIiK13ZUivH//p4+DB6Gg4Oq9jmM3tOjc2U6D6NTpk0fr1nZ5tGu99aZ9rSjIJ06cwOPxcPr0afr27cvUqVOJjIy8Cz+siEjVqCSLiNQWRUW2+H788aePyg/GhYbawtutG4wZA126XD06drSrRdym0tJSli1bxrp164iKimLu3Lncc889d/6ziYj4mUqyiEhNU1AAe/bYY/fuq6+HDtlR4yvatYPu3eHhh+Gee2wp7tbNFuGwML/HKioq4tVXX+X8+fMMHjyYSZMmEXEHhVtEpDqpJIuIBKvCQti7l4Ep93Om6BxdL4bb3eGKi6/e4zgQUw96REL9GKhf3x716lVshlEG7LLHaeyxxr8xfT4fm05sItQXCo3hscceo2PHjv79JiIifqaSLCIS6IyB48dh+3bYtu3q6/794PNdXTmiOBIaNbIlODLSHhERrq4RfPnyZc6dO0eoL5Tm9Zvz9NNPE1YNo9QiIv6mkiwiEkjKyuz0iC1brh7btkFe3tV7OnWC/v3hoYegb1+6nnkJ6tVj+YIV7uW+xqVLl0hJSWHn0Z20aNGCxMRE2rZt63YsEZEqU0kWEXFLaSns2gWZmfbYvNmOEl+ZLhERAf36wezZthT36wd9+9od5Cp76zd3P/sNGGPYuXMnixYtori4mHHjxjFq1ChCQ0PdjiYicktUkkVE7gafD/buhfXrr5bibduuFuJGjWDQIHjmGRg40B7du0Od4Plj+sKFCyxcuJB9+/bRtm1bEhMTadGihduxRERuS/D86SsiEkyys20h3rDhajG+st5wgwYweDA8+6x9jY21y6sF6SYaxhg2bdrE4sWLMcYwdepUhg4dSsj11ksWEQkSKskiIneqtBS2boW1ayEjw74eO2av1aljp0o88ggMHQrDhtnl1mpIgczNzcXr9XL06FE6depEQkICjRs3djuWiMgdU0kWEblVeXm2DK9ebY/MTLscG0DbtjBiBHzlKzB8OAwYYJdbq2F8Ph/r1q1j2bJlhIaGkpCQwMCBA7WltIjUGCrJIiI3c+IErFx5tRTv2GGXZatTx84d/uIX4d577dG+vdtpq92pU6fweDxkZ2fTvXt3ZsyYQYMGDdyOJSLiVyrJIiLXOnYMVqy4ehw4YM9HRdlR4pkzYdQoO30iMtLdrHdRWVkZq1atYvXq1URERDBz5kx69eql0WMRqZFUkkVETpyAZctg6VL7euSIPd+4MYwZA//2b/a1f/+gWm3Cn7KysvB4PJw5c4Z+/foxdepU6tev73YsEZFqUzv/tBeR2i0392opXroUPv7Ynm/cGMaNg+eeg7Fj7ZrENeQBu9tVUlLCsmXLWLduHQ0bNmTevHl069bN7VgiItVOJVlEar7iYlizBhYvtsfmzXZOcVSUHSH+whdgwgQ7UlzLS3Flhw4dwuv1kpeXR2xsLJMmTaJu3bpuxxIRuStUkkWk5jHGbu2ckgJpafahu8JCO1Vi+HD4/vdh0iQYMgTCwtxOG3CKiopIS0tjy5YtNGnShMcff5wOHTq4HUtE5K5SSRaRmuHCBUhPt8U4JeXqOsXdu8PnPw+TJ9spFNdu6SyfsHfvXhYuXMilS5cYOXIkY8eOJUz/kBCRWkglWUSCkzGwezcsXGiPjAwoK7O72U2aBN/9LkydChoBrZKLFy+SkpLCrl27aNmyJXPnzqVNmzZuxxIRcY1KsogEj8JC+8DdlWJ89Kg9378/fO1rEBdn1yrWyGeVGWPYvn07qamplJSUMGHCBEaMGEFoaKjb0UREXKWSLCKB7dQpSEoCj8c+dFdYCPXr29Hi73wHpk+Hdu3cThmU8vPzSUpK4sCBA7Rr147ExESaN2/udiwRkYCgkiwigeXKNAqPxx7r19tzMTHwxBOQkGDnFkdEuJ00aBljyMzMZMmSJRhjmDZtGkOGDCFEK3uIiPwflWQRcZ/PZ8vwBx/Y48oOd7Gx8IMfQGIi9OsH2tntjuXm5uLxeDh27BidO3cmISGB6Ohot2OJiAQclWQRcUdpKSxfDu+/Dx99BDk5dom2CRPgq1+1xbhtW7dT1hg+n4+MjAyWL19OWFgY9913H/3799eW0iIiN6CSLCJ3T3GxXabtvffgww/h/Hk7vzguDh54AGbMAI1q+t3JkyfxeDzk5OTQs2dP4uLiaNCggduxREQCmkqyiFSvoiK7ocd779k5xvn50KiRHSl+8EGYMgXq1XM7ZY1UVlbGypUrWbNmDfXq1WPWrFn06tXL7VgiIkFBJVlE/K+kxBbjv//dTqUoKIDGjeFzn4OZM2HiRND2xtXq+PHjeDwezp49S//+/Zk6dSr19I8REZEqq1JJdhznTSAeOG2M6XOd6w7wMjAduAw8bozZXHHtMeB7Fbf+yBjzJ38EF5EAU1oKS5faYvzBB5CXZ6dOzJoFDz1k5xpr/eJqV1JSQnp6Ohs2bKBRo0Y8/PDDdO3a1e1YIiJBp6ojyW8Bvwb+fIPrcUC3imMY8CowzHGcJsALQCxggE2O43iMMefvJLSIBAifD9auhb/9Df7xDzh71u54d//9MHu23Qo6PNztlLXGwYMH8Xq95OfnM2TIECZOnEhdjdiLiNyWKpVkY8xKx3E6fsYt9wF/NsYYYJ3jONGO47QGxgGLjTHnABzHWQxMA965k9Ai4rIdO2wxfucdu+tdvXp2/eK5c2HaNK1hfJf5fD7OnTvH22+/TdOmTVmwYAExMTFuxxIRCWr+mpPcFjhe6X1Wxbkbnf8Ux3GeAp4C9Ie7SCDKyoK//hXefht27oTQUPvQ3Y9+BPfdZ0eQ5a7bs2cPJ06cwOfzMWrUKMaOHUudOnrcRETkTvnrT9LrLbRpPuP8p08a8zrwOkBsbOx17xGRu+zCBbuO8V/+AsuW2Z3v7r0Xfv1rO89YWxi75uLFiyQnJ7Nnzx5Cw0Jp2bIlEydOdDuWiEiN4a+SnAW0r/S+HZBdcX7cNeeX++l7ikh1KC+HJUvgT3+yaxkXFkKXLvD88zB/PughMFcZY9i2bRupqamUlpYyceJElu5fqk1BRET8zF8l2QN82XGcd7EP7uUbY3Icx0kFfuI4TuOK+6YA3/bT9xQRf9qzxxbjv/wFsrPtkm2PPQaPPGJHj1XCXJeXl0dSUhIHDx4kJiaGhIQEmjVrhnNA/92IiPhbVZeAewc7ItzMcZws7IoVYQDGmN8Bydjl3w5gl4BbUHHtnOM4PwQ2VnypF688xCciASAvD959F956C9avt/OM4+Lg5Zftg3haGSEgGGPYsGED6enpOI5DXFwcQ4YM0eixiEg1qurqFnNvct0Az9zg2pvAm7ceTUSqhc8Hy5fDm2/Cv/5ld8Tr0wd+8Qt4+GFo1crthFLJ2bNn8Xg8HD9+nC5duhAfH0+0tu4WEal2egRapLY4ftyOGP/xj3D4sN0aesECeOIJGDxY0ykCTHl5ORkZGaxYsYLw8HDuv/9++vXrp9FjEZG7RCVZpCYrLYWkJPj97yElxa5OMWGCXbbtgQfs+sYScHJycvB4PJw8eZJevXoRFxdHVFSU27FERGoVlWSRmujgQfjDH+yo8cmT0KYNfPe7dtS4Uye308kNlJaWsmLFCjIyMoiMjOShhx6iZ8+ebscSEamVVJJFaoqSEvjoI3jtNUhPh5AQmDEDvvAF+zCeNpgIaMeOHcPj8ZCbm8uAAQOYMmUK9TTSLyLiGv2tKRLsjhyx0yn+8Ac4dQpiYuDFF+1843bt3E4nN1FcXEx6ejobN24kOjqaRx55hM6dO7sdS0Sk1lNJFglG5eWQnAy/+x0sWmQfupsxA770JZg61S7lJgHvwIEDJCUlkZ+fz7Bhw5gwYQLh4eFuxxIREVSSRYLL6dPwxht2SsWxY9C6NXzve/D5z9sRZAkKly9fJi0tjW3bttGsWTOeeOIJ2rdvf/NPFBGRu0YlWSTQGQNr18Jvfwv//KedezxxIvzyl5CYCGFhbieUKjLGsHv3bhYtWkRhYSFjxoxh9OjR1NF8cRGRgKM/mUUC1eXL8Le/wW9+A1u3QsOGdjrF009Djx5up5NbVFBQQHJyMnv37qV169bMnz+fVtq4RUQkYKkkiwSao0ftqPEbb8C5c9C3r517/PDDoLVyg44xhq1bt5Kamkp5eTmTJk3i3nvvJSQkxO1oIiLyGVSSRQKBMXar6FdeAY/HPoj3wAPw7LMwerR2wwtS58+fJykpiUOHDtGhQwcSEhJo2rSp27FERKQKVJJF3FRYCH/9K7z8MuzcCU2bwje/aadU6EGuoOXz+diwYQNLly7FcRxmzJjB4MGDtaW0iEgQUUkWcUN2tp1S8bvfQW4uDBgAb74Jc+dCRITb6eQOnDlzBo/HQ1ZWFt26dWPGjBk0atTI7VgiInKLVJJF7qaNG+F//gf+8Q+71vF998FXvgJjxmhKRZArLy9n9erVrFq1ivDwcB544AH69u2r0WMRkSClkixS3crL7TzjX/4SVq+GBg3sXOMvfxm0s1qNkJ2djcfj4dSpU/Tp04dp06YRGRnpdiwREbkDKski1eXiRXjrLTtyfPAgdOwIv/oVPPGEXc5Ngl5paSnLly9n7dq1REVFMWfOHLp37+52LBER8QOVZBF/y862q1S89hrk5cHw4fDTn8L994M2jagxjhw5gtfr5dy5cwwaNIjJkycTofnkIiI1hv7GFvGXXbvgF7+wq1WUl8PnPgdf/Srce6/bycSPiouLWbx4MZs2baJx48Y8+uijdOrUye1YIiLiZyrJInfiyvrGL70EixZBvXrwxS/Cc89pvnENtG/fPhYuXEhBQQHDhw9nwoQJhGlbcBGRGkklWeR2lJfDv/5ly3FmJjRvDi++CP/2b3atY6lRLl++TEpKCjt27KB58+bMmjWLdu3auR1LRESqkUqyyK0oKoI//clOqzhwALp1s2sdP/qoHUWWGsUYw65du1i0aBFFRUWMHTuW0aNHExoa6nY0ERGpZirJIlWRlwevvmp3xjt1CoYMgffesw/jqTDVSBcuXCA5OZmPP/6YNm3akJiYSMuWLd2OJSIid4lKsshnOXnSLtv26qtQUABTp9pto8eN0+YfNZQxhs2bN7N48WLKy8uZPHkyw4cPJyQkxO1oIiJyF6kki1zPoUN2vvEf/wilpfDQQ7YcDxjgdjKpRufOncPr9XLkyBE6duxIQkICTZo0cTuWiIi4QCVZpLKdO+2axu++a6dRPP44fP3r0LWr28mkGvl8PtavX8/SpUsJDQ0lPj6eQYMGaUtpEZFaTCVZBGDjRvjxj+GjjyAyEr7yFbvGcZs2bieTanb69Gk8Hg8nTpzgnnvuYcaMGTTUjogiIrWeSrLUbitX2nKclgaNG8MLL8Czz2oZt1qgvLycVatWsWrVKiIiInjwwQfp3bu3Ro9FRARQSZbayBhbin/0I1i9Glq0gJ/9DJ5+Gho0cDud3AUnTpzA4/Fw+vRp+vbty7Rp06hfv77bsUREJICoJEvtYQwkJdlNPzIzoV07eOUVePJJUEGqFUpLS1m2bBnr1q0jKiqKuXPncs8997gdS0REApBKstR8Ph98+CH88IewdSt06gSvvQaPPQZ167qdTu6Sw4cP4/V6OX/+PIMHD2bSpElERES4HUtERAKUSrLUXFe2jv7hD+2qFV272iXdHn4YwsLcTid3SVFREYsXL2bz5s00adKExx57jI4dO7odS0REApxKstQ85eXwz3/aaRV79kCPHvD22zB7NtTRr3xt8vHHH7Nw4UIuXrzIiBEjGDduHGH6B5KIiFSBGoPUHNeW49697XrHM2dq6+ha5tKlS6SkpLBz505atGjBnDlzaKPl/ERE5BaoJEvwu145/sc/4MEHQVsJ1yrGGHbu3MmiRYsoLi5m3LhxjBo1ilD9I0lERG6RSrIEL5/PluMf/EDlWMjPz2fhwoXs37+fdu3akZCQQIsWLdyOJSIiQUolWYKPzwcffGA3/ti1C3r1gr//3U6rUDmudYwxbNq0icWLF2OMYerUqQwdOpQQ/S6IiMgdUEmW4GGM3Tb6+9+Hbduge3d45x2YNUtzjmup3NxcvF4vR48epVOnTiQkJNC4cWO3Y4mISA2gkiyBzxhITobnn4fNm6FbN7taxZw5Kse1lM/nY+3atSxfvpzQ0FASExMZMGCAtpQWERG/UUmWwGUMLF0K3/serFtnNwF56y27zrGWcqu1Tp06hcfjITs7m+7duzNjxgwaaDtxERHxMzUNCUyrV8P/+3+wfLndPvq112DBAm0CUouVlZWxatUqVq9eTb169Zg5cya9evXS6LGIiFSLKpVkx3GmAS8DocAbxpifXnP9V8D4irf1gRbGmOiKa+XAjoprx4wxif4ILjVUZqYtxykp0LIlvPIKfOELoO2Da7WsrCw8Hg9nzpyhX79+TJ06lfr167sdS0REarCblmTHcUKB3wCTgSxgo+M4HmPM7iv3GGOeq3T/s8DASl+i0BgzwH+RpUbatcvOOX7/fWjSBH7+c3jmGVARqtVKSkpYunQp69evp2HDhsybN49u3bq5HUtERGqBqowkDwUOGGMOATiO8y5wH7D7BvfPBV7wTzyp8Q4dsqtVvP02REXZj597Dho2dDuZuOzQoUN4vV7y8vIYMmQIEydOpG7dum7HEhGRWqIqJbktcLzS+yxg2PVudBynA9AJWFrpdITjOJlAGfBTY8yHt5lVapLsbPjhD+GNN+xDeF/7Gnzzm9C0qdvJxGVFRUWkpqaydetWmjRpwuOPP06HDh3cjiUiIrVMVUry9Z6KMTe4dw7wnjGmvNK5GGNMtuM4nYGljuPsMMYc/NQ3cZyngKcAYmJiqhBLgtK5c/Czn9m5xuXl8NRT8N3vQps2bieTALB3714WLlzIpUuXGDlyJGPHjiVMD2uKiIgLqlKSs4D2ld63A7JvcO8c4JnKJ4wx2RWvhxzHWY6dr/ypkmyMeR14HSA2NvZGJVyC1aVL8PLLdq7xhQswf77dTrpTJ7eTSQC4ePEiixYtYvfu3bRs2ZK5c+fSRv9wEhERF1WlJG8EujmO0wk4gS3C8669yXGc7kBjYG2lc42By8aYYsdxmgEjgZ/7I7gEiZIS+P3v7dSKU6cgMRF+9CPo29ftZBIAjDFs376d1NRUSkpKmDBhAiNGjCBUm8SIiIjLblqSjTFljuN8GUjFLgH3pjFml+M4LwKZxhhPxa1zgXeNMZVHgXsCrzmO4wNCsHOSb/TAn9QkPh+8+65dzu3QIRgzxq5cMWKE28kkQOTn55OUlMSBAwdo3749iYmJNGvWzO1YIiIiQBXXSTbGJAPJ15x7/pr337/O52UAGjKsTYyB1FT49rdh61bo399uKT1tGmjTB8GOHm/cuJH09HSMMUybNo2hQ4dqUxAREQko2nFP/GfDBvjWt2DZMjvX+K9/hTlzICTE7WQSIM6ePYvX6+XYsWN07tyZhIQEoqOj3Y4lIiLyKSrJcuf27bMrVLz3HjRvbleu+OIXITzc7WQSIHw+HxkZGSxfvpywsDDuu+8++vfvr9FjEREJWCrJcvtOnbIrVLz+ut02+oUX4D//Exo0cDuZBJCTJ0/i8XjIycmhZ8+eTJ8+naioKLdjiYiIfCaVZLl1BQXw3/8Nv/gFFBfbUePnn4eWLd1OJgGkrKyMFStWsGbNGurXr8+sWbPo1auX27FERESqRCVZqq601C7n9oMfwOnTMGsW/PjH0K2b28kkwBw7dgyv18vZs2cZMGAAU6ZMoV69em7HEhERqTKVZLk5Y+CDD+xDefv3w9ix4PXC0KFuJ5MAU1JSQnp6Ohs2bKBRo0bMnz+fLl26uB1LRETklqkky2fLyICvf92+9uxpy/GMGVrOTT7l4MGDeL1e8vPzGTp0KBMnTiRcD2+KiEiQUkmW69u/344cv/8+tGplH85bsADq6FdGPqmwsJC0tDS2bt1K06ZNWbBgATExMW7HEhERuSNqPPJJZ8/Ciy/Cq6/aFStefBG++lWIjHQ7mQSg3bt3k5yczOXLlxk1ahRjx46ljv4hJSIiNYD+NhOrqAhefhl+8hO4dAm+8AX4/ve1YoVc18WLF0lOTmbPnj20atWK+fPn06pVK7djiYiI+I1Kcm3n88E778B3vgPHjkF8PPzsZ6CluuQ6jDFs27aN1NRUSktLmThxIiNGjCBEuyqKiEgNo5Jcm61caTf/yMyEQYPgrbdg/Hi3U0mAysvLw+v1cujQIWJiYkhISKBZs2ZuxxIREakWKsm10f798I1vwIcfQvv28Je/wLx5oNFAuQ5jDBs2bCA9PR3HcZg+fTqxsbHaUlpERGo0leTa5Nw5+yDeb35jH8r78Y/huedAmzzIDZw5cwav18vx48fp2rUrM2bMIDo62u1YIiIi1U4luTYoKbHF+MUX4cIF+Pzn7a55etBKbqC8vJyMjAxWrFhBeHg4999/P/369dPosYiI1BoqyTWZMXZKxde/DgcPwpQp8N//DX36uJ1MAlhOTg4ej4eTJ0/Su3dvpk2bRlRUlNuxRERE7iqV5Jpqyxa7vvHy5XaliuRkiItzO5UEsNLSUlasWEFGRgaRkZHMnj2bHj16uB1LRETEFSrJNU12Nnz3u/CnP0HTpnaaxVNPaac8+UxHjx7F6/WSm5vLwIEDmTx5MvU0V11ERGoxNaeaorDQTqX46U+htBS+9jW79rEespLPUFxczJIlS8jMzCQ6OppHHnmEzp07ux1LRETEdSrJwc4Y+Pvf7ZJux4/D5z4HP/85dOnidjIJcPv37ycpKYkLFy4wbNgwJkyYQHh4uNuxREREAoJKcjDbsMEu4ZaRAQMH2vWOx451O5UEuMuXL5Oamsr27dtp3rw5TzzxBO3bt3c7loiISEBRSQ5GJ07At79tS3HLlvCHP8Bjj0FoqNvJJIAZY9i9ezeLFi2isLCQMWPGMHr0aOpovrqIiMin6G/HYFJYCL/4hZ13XFYG3/qWnXfcoIHbySTA/f/27jw6yvu+9/j7i9jBiM0YDMjsxoANxgIMGBsMCBBCJN6QbWy8JD7Nnpv2ts1tj5s6bY6TnJOtSVu7iXvSnps4dW5az2hBLGKxWQUYg4GAMfu+iH3V8r1/zJAzmQhrZEt6ZjSf1zk6zLOh73z5zcNHj57lwoULFBcX8/vf/55evXrxzDPPcNtttwVdloiISNJSSE4F7vDWW5Hzjvfvh0cfhe9/H/r3D7oySXLuznvvvceiRYuorq5m2rRpjB8/nhZ6BLmIiMjHUkhOdps2wde/Du+8A/fcA8uWweTJQVclKeDMmTOEw2H27t3LHXfcwZw5c+jWrVvQZYmIiKQEheRkdeJE5H7Hv/hF5H7Hr70GL76o846lTjU1Naxfv56ysjLMjNmzZ3PffffpkdIiIiL1oJCcbK5fh3/6J3jlFbh8OXIU+eWXdb9jScjJkycJhUIcOnSIwYMHM3v2bDIzM4MuS0REJOUoJCeTkpLILd127ow8QvqHP4Q77wy6KkkB1dXVvPvuu7zzzju0bt2aRx55hBEjRujosYiIyCekkJwMdu2KhOPiYhgyBIqKIDc36KokRRw+fJhQKMSJEycYMWIEM2fOpEOHDkGXJSIiktIUkoN0/jx8+9vw4x9D27aR27t95Sugp55JAiorK1m+fDlr1qyhY8eOFBQUcKd+8yAiItIgFJKDUFMTeRDIX/0VHD8OL7wA3/lO5MEgIgnYt28f4XCYiooKRo8ezfTp02nbtm3QZYmIiDQbCslNrbw8crR43Tq4/34Ih2HMmKCrkhRx7do1Fi9ezMaNG+nSpQvPPvss/XW/bBERkQankNxUjh+PPEr63/8devaEX/4S5s8HPdRBErRr1y6Kioq4cOEC48ePZ8qUKbRq1SroskRERJolheTGVlkJP/0pfOtbkcdK/+Vfwt/+rR4lXpHPJAAAGolJREFULQm7dOkSpaWlbN26lR49evDEE0/Qu3fvoMsSERFp1hSSG9OSJfDVr8KOHTBzJvzoR7qlmyTM3dm2bRslJSVcvXqVhx56iEmTJpGhB8qIiIg0OoXkxrBvH/z5n8PvfgcDBkAoBHl5oHvWSoLOnz9PUVERu3btonfv3uTn59OjR4+gyxIREUkbCskN6coV+N734NVXI+ca/+M/wje+Ebm9m0gC3J1NmzaxePFiqqurycnJYdy4cbTQuesiIiJNSiG5IbhHjhZ//euRo8hPPBG553HfvkFXJimkoqKCcDjMvn376NevH3PmzKFr165BlyUiIpKWFJI/rZ074Wtfg9JSGD4cyspgypSgq5IUUlNTw7p16ygrKyMjI4O8vDxGjx6tR0qLiIgESCH5k7p4Ef7hH+AHP4B27SIX5X3xi6Bbckk9nDhxglAoxOHDhxkyZAizZ8+mU6dOQZclIiKS9hI60dHMZprZTjPbbWZ/Xcvy58zspJltjn59LmbZAjP7MPq1oCGLD4Q7/OY3MHQofPe7kXsd79oVOZqsgCwJqq6uZvny5bz22mucOXOGRx99lIKCAgVkERGRJFHnkWQzywB+BkwHDgHlZhZy9+1xq/7G3b8ct21X4O+AbMCBjdFtzzRI9U1t27bI0/KWLYPRo+Gtt2D8+KCrkhRz+PBh3n77bU6ePMndd9/NzJkzad++fdBliYiISIxETrcYC+x29z0AZvYmMBeID8m1mQEsdveK6LaLgZnArz9ZuQE5fx7+/u/hJz+JPATkX/4FPv950P1qpR4qKyspKytj3bp13HLLLTz55JMMGTIk6LJERESkFomE5N7AwZjpQ8C4WtZ71MweBHYB/8vdD95k29R5VJg7/OpX8Bd/EXms9Oc+B9/5DnTvHnRlkmL27t1LOBzmzJkz3HfffUyfPp02bdoEXZaIiIjcRCIhubZL7D1uOgz82t2vmdmfAb8EHk5w28g3MXsJeAkgKysrgbIa2QcfwJe+BCtXwpgx8PbbMHZs0FVJirl69SqLFy9m06ZNdO3alQULFtCvX7+gyxIREZE6JBKSDwGxN/ztAxyJXcHdT8dM/hvw3ZhtJ8dtu7y2b+LurwOvA2RnZ9capJvE+fPwrW9FTq3IzITXX4cXX4w8HESkHnbu3ElRUREXL15kwoQJTJ48mVa6uFNERCQlJBKSy4HBZtYfOAwUAE/FrmBmvdz9aHQyH9gRfV0KfMfMukSnc4BvfuqqG0P8qRWf/3zk1Ipu3YKuTFLMpUuXKCkpYdu2bfTo0YOCggJuv/32oMsSERGReqgzJLt7lZl9mUjgzQDecPdtZvYKsMHdQ8BXzSwfqAIqgOei21aY2beJBG2AV25cxJdU3CEvD4qLITs78vS8MWOCrkpSjLuzdetWFi5cyPXr15kyZQoTJ04kQxd4ioiIpJyEHibi7sVAcdy8l2Nef5ObHCF29zeANz5FjY3PLBKS586NnFqhUCP1dO7cOYqKivjwww/p06cP+fn53HrrrUGXJSIiIp+Qnrh3wxe+EHQFkoLcnY0bN7J48WLcnRkzZjB27Fha6Bx2ERGRlKaQLPIJnT59mnA4zP79+xkwYAB5eXl06dKl7g1FREQk6Skki9RTTU0Na9asYfny5WRkZJCfn8+oUaMwq+2OhyIiIpKKFJJF6uHYsWOEQiGOHj3K0KFDyc3N5ZZbbgm6LBEREWlgCskiCaiqqmLlypWsWrWKdu3a8dhjjzFs2DAdPRYREWmmFJJF6nDw4EFCoRCnTp1i5MiR5OTk0L59+6DLEhERkUakkCxyE9evX6esrIx169aRmZnJ008/zaBBg4IuS0RERJqAQrJILfbs2UM4HObs2bOMGTOGqVOn0qZNm6DLEhERkSaikCwS48qVKyxatIjNmzfTrVs3nnvuOe64446gyxIREZEmppAsErVjxw6Ki4u5dOkSEydOZPLkybRsqY+IiIhIOlICkLR38eJFSkpK2L59Oz179uSpp56iV69eQZclIiIiAVJIlrTl7mzZsoWFCxdSWVnJww8/zIQJE8jIyAi6NBEREQmYQrKkpXPnzlFYWMju3bvp27cv+fn5dO/ePeiyREREJEkoJEtacXfKy8tZunQp7s6sWbMYM2aMHgoiIiIif0QhWdLGqVOnCIfDHDhwgIEDB5KXl0fnzp2DLktERESSkEKyNHvV1dWsWbOG5cuX06pVK+bOncvIkSN19FhERERuSiFZmrWjR48SCoU4duwYd911F7m5uXTs2DHoskRERCTJKSRLs1RVVcWKFStYtWoV7du35/HHH2fYsGFBlyUiIiIpQiFZmp0DBw4QCoU4ffo0o0aNIicnh3bt2gVdloiIiKQQhWRpNq5fv87SpUtZv349mZmZzJ8/n4EDBwZdloiIiKQghWRpFnbv3k1hYSHnzp1j7NixTJ06ldatWwddloiIiKQohWRJaVeuXKG0tJT333+f7t278/zzz5OVlRV0WSIiIpLiFJIlZW3fvp3i4mIuX77MpEmTePDBB2nZUkNaREREPj0lCkk5Fy5coKSkhB07dtCrVy/mz59Pz549gy5LREREmhGFZEkZ7s77779PaWkplZWVTJ06lQkTJtCiRYugSxMREZFmRiFZUsLZs2cJh8Ps2bOHrKws8vPz6datW9BliYiISDOlkCxJraamhvLycpYuXYqZkZubS3Z2th4pLSIiIo1KIVmS1smTJwmHwxw8eJBBgwaRl5dHZmZm0GWJiIhIGlBIlqRTXV3NqlWrWLlyJa1bt+Yzn/kM99xzj44ei4iISJNRSJakcvToUd5++22OHz/O8OHDmTlzJh07dgy6LBEREUkzCsmSFCorK1mxYgWrV6+mQ4cOzJs3j6FDhwZdloiIiKQphWQJ3P79+wmFQlRUVHDvvfeSk5ND27Ztgy5LRERE0phCsgTm2rVrLFmyhA0bNtC5c2eeeeYZBgwYEHRZIiIiIgrJEowPP/yQwsJCzp8/z7hx43j44Ydp3bp10GWJiIiIAArJ0sQuX75MaWkpW7Zs4dZbb+XFF1+kT58+QZclIiIi8kcUkqVJuDvbt2+nuLiYq1ev8uCDDzJp0iRattQQFBERkeSjhCKN7sKFCxQVFbFz505uv/128vPzue2224IuS0REROSmFJKl0bg77733HosWLaK6uprp06dz//3306JFi6BLExEREflYCsnSKM6cOUM4HGbv3r3ccccd5Ofn07Vr16DLEhEREUmIQrI0qJqaGtavX09ZWRlmxuzZs7nvvvv0SGkRERFJKQrJ0mBOnDhBKBTi8OHDDB48mLy8PDp16hR0WSIiIiL1llBINrOZwI+BDODn7v5q3PJvAJ8DqoCTwAvuvj+6rBrYGl31gLvnN1DtkiSqq6t59913WblyJW3atOGRRx5hxIgROnosIiIiKavOkGxmGcDPgOnAIaDczELuvj1mtfeAbHe/bGZfAL4HzIsuu+Luoxq4bkkShw8fJhQKceLECUaMGMHMmTPp0KFD0GWJiIiIfCqJHEkeC+x29z0AZvYmMBf4Q0h292Ux668F5jdkkZJ8KisrWbZsGWvXrqVjx44UFBRw5513Bl2WiIiISINIJCT3Bg7GTB8Cxn3M+i8CJTHTbc1sA5FTMV519/+pd5WSVPbt20c4HKaiooLRo0czffp02rZtG3RZIiIiIg0mkZBc24mlXuuKZvOBbOChmNlZ7n7EzAYAZWa21d0/qmXbl4CXALKyshIoS5ra1atXWbJkCRs3bqRLly48++yz9O/fP+iyRERERBpcIiH5ENA3ZroPcCR+JTObBvwN8JC7X7sx392PRP/cY2bLgXuBPwnJ7v468DpAdnZ2rSFcgrNr1y4KCwu5ePEi48ePZ8qUKbRq1SroskREREQaRSIhuRwYbGb9gcNAAfBU7Apmdi/wGjDT3U/EzO8CXHb3a2bWHZhI5KI+SRGXLl2itLSUrVu30qNHD+bNm0fv3r2DLktERESkUdUZkt29ysy+DJQSuQXcG+6+zcxeATa4ewj4PtAReCt6268bt3q7C3jNzGqAFkTOSd5e6zeSpOLufPDBByxcuJCrV68yefJkHnjgATIyMoIuTURERKTRJXSfZHcvBorj5r0c83raTbZbDdz9aQqUpnf+/HmKiorYtWsXvXv3Jj8/nx49egRdloiIiEiT0RP35A/cnU2bNrF48WKqq6vJyclh3LhxtGjRIujSRERERJqUQrIAUFFRQTgcZt++ffTr1485c+bQtWvXoMsSERERCYRCcpqrqalh7dq1LFu2jIyMDObMmcO9996rR0qLiIhIWlNITmPHjx8nFApx5MgR7rzzTnJzc+nUqVPQZYmIiIgETiE5DVVVVfHOO+/w7rvv0rZtWx599FGGDx+uo8ciIiIiUQrJaebQoUOEQiFOnjzJPffcw4wZM2jfvn3QZYmIiIgkFYXkNHH9+nWWLVvG2rVr6dSpE08++SRDhgwJuiwRERGRpKSQnAb27t1LOBzmzJkzZGdnM23aNNq0aRN0WSIiIiJJSyG5Gbt69SqLFi3ivffeo2vXrixYsIB+/foFXZaIiIhI0lNIbqZ27txJUVERFy9eZMKECUyePJlWrVoFXZaIiIhISlBIbmYuXbpESUkJ27Zt47bbbqOgoIDbb7896LJEREREUopCcjPh7mzdupWFCxdy/fp1pkyZwsSJE8nIyAi6NBEREZGUo5DcDJw7d46ioiI+/PBD+vTpQ35+PrfeemvQZYmIiIikLIXkFObubNiwgSVLluDuzJgxg7Fjx9KiRYugSxMRERFJaQrJKer06dOEw2H279/PgAEDyMvLo0uXLkGXJSIiItIsKCSnmJqaGtasWcPy5ctp2bIl+fn5jBo1So+UFhEREWlACskp5NixY4RCIY4ePcrQoUPJzc3llltuCbosERERkWZHITkFVFVVsXLlSlatWkW7du14/PHHueuuu3T0WERERKSRKCQnuYMHDxIKhTh16hQjR44kJyeH9u3bB12WiIiISLOmkJykrl+/TllZGevWrSMzM5Onn36aQYMGBV2WiIiISFpQSE5CH330EYWFhZw9e5YxY8YwdepU2rRpE3RZIiIiImlDITmJXLlyhUWLFrF582a6devG888/T1ZWVtBliYiIiKQdheQksWPHDoqLi7l06RIPPPAADz30EC1b6p9HREREJAhKYQG7ePEiJSUlbN++nZ49e/LUU0/Rq1evoMsSERERSWsKyQFxd7Zs2cLChQuprKzk4YcfZsKECWRkZARdmoiIiEjaU0gOwNmzZyksLOSjjz6ib9++5Ofn071796DLEhEREZEoheQm5O6Ul5ezdOlS3J1Zs2YxZswYPRREREREJMkoJDeRU6dOEQ6HOXDgAAMHDiQvL4/OnTsHXZaIiIiI1EIhuZFVV1ezevVqVqxYQatWrZg7dy4jR47U0WMRERGRJKaQ3IiOHj1KKBTi2LFjDBs2jFmzZtGxY8egyxIRERGROigkN4KqqipWrFjBqlWraN++PU888QR33XVX0GWJiIiISIIUkhvYgQMHCIVCnD59mlGjRpGTk0O7du2CLktERERE6kEhuYFcu3aNpUuXUl5eTufOnZk/fz4DBw4MuiwRERER+QQUkhvA7t27KSws5Ny5c4wdO5apU6fSunXroMsSERERkU9IIflTuHLlCqWlpbz//vt0796dF154gb59+wZdloiIiIh8SgrJn9D27dspLi7mypUrTJo0iQcffJCWLdVOERERkeZAqa6eLly4QElJCTt27KBXr17Mnz+fnj17Bl2WiIiIiDQgheQEuTubN29m0aJFVFVVMW3aNMaPH0+LFi2CLk1EREREGphCcgLOnDlDYWEhe/bsISsri/z8fLp16xZ0WSIiIiLSSBSSP0ZNTQ3l5eUsXboUMyM3N5fs7Gw9UlpERESkmUsoJJvZTODHQAbwc3d/NW55G+A/gPuA08A8d98XXfZN4EWgGviqu5c2WPWN6OTJk4RCIQ4dOsSgQYPIy8sjMzMz6LJEREREpAnUGZLNLAP4GTAdOASUm1nI3bfHrPYicMbdB5lZAfBdYJ6ZDQMKgOHA7cASMxvi7tUN/UYaSnV1NatWrWLlypW0bt2az372s9x99906eiwiIiKSRhI5kjwW2O3uewDM7E1gLhAbkucC34q+/i3wU4ukyrnAm+5+DdhrZrujf9+ahim/YR05coRQKMTx48cZPnw4s2bNokOHDkGXJSIiIiJNLJGQ3Bs4GDN9CBh3s3XcvcrMzgHdovPXxm3b+xNX20jcnaVLl7J69Wo6dOjAvHnzGDp0aNBliYiIiEhAEgnJtZ1n4Amuk8i2kb/A7CXgJYCsrKwEymo4Zsa1a9cYNWoUOTk5tG3btkm/v4jIpzGq56igSxARaXYSCcmHgNhnLfcBjtxknUNm1hLIBCoS3BYAd38deB0gOzu71iDdmHJzc3XesYikpB/N/FHQJYiINDuJPAmjHBhsZv3NrDWRC/FCceuEgAXR148BZe7u0fkFZtbGzPoDg4H1DVN6w1JAFhEREZEb6jySHD3H+MtAKZFbwL3h7tvM7BVgg7uHgF8A/xm9MK+CSJAmut5/EbnIrwr4UjLf2UJEREREBMAiB3yTS3Z2tm/YsCHoMkRERESkGTOzje6eXduyRE63EBERERFJKwrJIiIiIiJxFJJFREREROIoJIuIiIiIxFFIFhERERGJo5AsIiIiIhJHIVlEREREJI5CsoiIiIhIHIVkEREREZE4CskiIiIiInEUkkVERERE4igki4iIiIjEUUgWEREREYmjkCwiIiIiEkchWUREREQkjrl70DX8CTM7CewP4Ft3B04F8H1TlfpVP+pX/ahf9aN+1Y/6VX/qWf2oX/UTVL/ucPdba1uQlCE5KGa2wd2zg64jVahf9aN+1Y/6VT/qV/2oX/WnntWP+lU/ydgvnW4hIiIiIhJHIVlEREREJI5C8h97PegCUoz6VT/qV/2oX/WjftWP+lV/6ln9qF/1k3T90jnJIiIiIiJxdCRZRERERCRO2oRkM5tpZjvNbLeZ/XUty9uY2W+iy9eZWb+YZd+Mzt9pZjOasu6gJNCvb5jZdjPbYmZLzeyOmGXVZrY5+hVq2sqDkUC/njOzkzF9+VzMsgVm9mH0a0HTVh6MBPr1w5he7TKzszHL0mp8mdkbZnbCzD64yXIzs59Ee7nFzEbHLEvHsVVXv56O9mmLma02s5Exy/aZ2dbo2NrQdFUHK4GeTTazczGfu5djln3sZ7k5SqBf/zumVx9E91ldo8vSaoyZWV8zW2ZmO8xsm5l9rZZ1kncf5u7N/gvIAD4CBgCtgfeBYXHrfBH41+jrAuA30dfDouu3AfpH/56MoN9TEvRrCtA++voLN/oVnb4Y9HtIwn49B/y0lm27Anuif3aJvu4S9HsKul9x638FeCNmOt3G14PAaOCDmyzPBUoAA+4H1qXr2EqwXxNu9AGYdaNf0el9QPeg30MS9mwyUFjL/Hp9lpvLV139ilt3DlAWM51WYwzoBYyOvr4F2FXL/49Juw9LlyPJY4Hd7r7H3a8DbwJz49aZC/wy+vq3wFQzs+j8N939mrvvBXZH/77mrM5+ufsyd78cnVwL9GniGpNJIuPrZmYAi929wt3PAIuBmY1UZ7Kob7+eBH7dJJUlIXdfCVR8zCpzgf/wiLVAZzPrRXqOrTr75e6ro/0A7buAhMbYzXyafV/Kqme/0n3/ddTdN0VfXwB2AL3jVkvafVi6hOTewMGY6UP86T/SH9Zx9yrgHNAtwW2bm/q+5xeJ/BR4Q1sz22Bma83sM41RYJJJtF+PRn+V9Fsz61vPbZuThN9z9DSe/kBZzOx0G191uVk/03Fs1Vf8vsuBRWa20cxeCqimZDXezN43sxIzGx6dpzH2McysPZFQ9/9iZqftGLPIaaz3AuviFiXtPqxlU36zAFkt8+Jv63GzdRLZtrlJ+D2b2XwgG3goZnaWux8xswFAmZltdfePGqHOZJFIv8LAr939mpn9GZHfWjyc4LbNTX3ecwHwW3evjpmXbuOrLtp3fQJmNoVISH4gZvbE6NjqASw2s99Hjxqmu01EHt170cxygf8BBqMxVpc5wCp3jz3qnJZjzMw6Evlh4evufj5+cS2bJMU+LF2OJB8C+sZM9wGO3GwdM2sJZBL5dUoi2zY3Cb1nM5sG/A2Q7+7Xbsx39yPRP/cAy4n85Nic1dkvdz8d06N/A+5LdNtmqD7vuYC4X1Wm4fiqy836mY5jKyFmdg/wc2Cuu5++MT9mbJ0A/pvmf2pdQtz9vLtfjL4uBlqZWXc0xurycfuvtBljZtaKSED+v+7+u1pWSdp9WLqE5HJgsJn1N7PWRAZu/FXxIeDGlZOPETnR3qPzCyxy94v+RH56Xt9EdQelzn6Z2b3Aa0QC8omY+V3MrE30dXdgIrC9ySoPRiL96hUzmU/kvCyAUiAn2rcuQE50XnOWyOcRM7uTyMUaa2LmpeP4qksIeDZ6hfj9wDl3P0p6jq06mVkW8DvgGXffFTO/g5ndcuM1kX7VeveCdGNmPaPX6GBmY4lkh9Mk+FlOR2aWSeQ3rG/HzEu7MRYdN78Adrj7D26yWtLuw9LidAt3rzKzLxNpbgaRK+W3mdkrwAZ3DxH5R/xPM9tN5AhyQXTbbWb2X0T+I64CvhT3q99mJ8F+fR/oCLwV3XcecPd84C7gNTOrIbIjfdXdm3WISbBfXzWzfCJjqILI3S5w9woz+zaR/2wAXon71Vyzk2C/IHLBy5vRH1ZvSLvxZWa/JnJ3ge5mdgj4O6AVgLv/K1BM5Orw3cBl4PnosrQbW5BQv14mcr3JP0f3XVXung3cBvx3dF5L4FfuvrDJ30AAEujZY8AXzKwKuAIURD+XtX6WA3gLTSqBfgF8Fljk7pdiNk3HMTYReAbYamabo/P+D5AFyb8P0xP3RERERETipMvpFiIiIiIiCVNIFhERERGJo5AsIiIiIhJHIVlEREREJI5CsoiIiIhIHIVkEREREZE4CskiIiIiInEUkkVERERE4vx/Z4v8xVOuOVUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "a,b,x0 = 0,2,1.0  # bounds and starting value\n",
    "# a,b,x0 = 0,2.5,2.75\n",
    "# a,b,x0 = 0,5,2.85\n",
    "xd = np.linspace(a,b,100)  # x grid\n",
    "plt.plot([a,b],[a,b],c='grey')  # plot the 45 degree line\n",
    "plt.plot(xd,F(xd),c='red')  # plot the function\n",
    "def plot_step(**kwargs):\n",
    "    if plot_step.counter < 10:\n",
    "        if plot_step.counter == 0:\n",
    "            x,f = kwargs['x0'],F(kwargs['x0'])\n",
    "            plt.plot([x,x],[0,f],c='green')  # initial vertical line\n",
    "            plt.plot([x,f],[f,f],c='green')\n",
    "        plot_step.counter += 1\n",
    "        x,f = kwargs['x'],F(kwargs['x'])\n",
    "        plt.plot([x,x],[x,f],c='green')  # vertical line\n",
    "        plt.plot([x,f],[f,f],c='green')  # horizontal line\n",
    "plot_step.counter = 0  # new public attribute\n",
    "x = solve_sa(F,x0,tol=1e-10,callback=plot_step)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Condition for convergence in scalar case\n",
    "\n",
    "$$\n",
    "|F'(x^\\star)| < 1\n",
    "$$\n",
    "\n",
    "- where $ x^\\star $ is the solution/root  \n",
    "- **stable solution**  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAHSCAYAAAAezFYoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd5iU1d3G8e+hKwoqYhdRY+yCumKJsaGIYIkRFCFqbNjQ2LuimMTeoqjB3hskgiACSywRRVksBAsWBCVGJYgaVPrz/nGW1xUXWGB2z5Tv57rmmt2Z2d0bE+WeZ3/nnJBlGZIkSZJ+VC91AEmSJCnfWJIlSZKkhViSJUmSpIVYkiVJkqSFWJIlSZKkhViSJUmSpIU0SB2gOquvvnrWunXr1DEkSZJUxMaOHfvfLMtaVvdcXpbk1q1bU1FRkTqGJEmSilgIYfKinnPcQpIkSVqIJVmSJElaiCVZkiRJWoglWZIkSVqIJVmSJElaiCVZkiRJWoglWZIkSVqIJVmSJElaiCVZkiRJWoglWZIkSVqIJVmSJElaiCVZkiRJWoglWZIkSVqIJVmSJElaiCVZkiRJWoglWZIkSVpIg9QBJEkFYPZsmDMH6teHJk1Sp5GkWmdJliT91NSpMHgwjBoFY8fC5MkwffqPzzdtCq1awdZbw847w377wS9/CSGkyyxJOWZJliTB3LkwaBDcfjv84x8wfz6sthqUlcEuu8Daa0PjxvFq8n//CxMnwpgx8MQTcMYZ8ItfwDHHxNuaa6b+00jScrMkS1Ipmz8/Ft1LL4UPPoD114eLLoKDD4a2bZd8dXjSJBg6NH6PCy+M3+d3v4PevaF167r4E0hSrXDhniSVqgkTYPfd4fDD41Xi/v3jFeI+fWDbbWs2PtG6NZx0Ejz3HLz3Hpx8Mjz6aBy/6NXrp2MaklRALMmSVGrmz4drr4U2bWD8eLj7bnjzTTjkEGiwHL9g3HRTuPlm+OijOHZx++2w+ebw+OOQZbnLL0l1wJIsSaXk229jGT73XOjUCd59Nxba+vVz9zPWXRfuuCPOLK+3HnTrBoceCl9/nbufIUm1zJIsSaXio49gxx3h6afhxhthwABYa63a+3nbbQejR8OVV8JTT8UZ51deqb2fJ0k5ZEmWpFIwfjzsuit8+SWUl8Ppp9fNlm0NGsD558NLL8Wft/vucO+9tf9zJWk5WZIlqdiNGRPLab168M9/wh571H2GHXeE11+POY45Bs4+G+bNq/scklRDlmRJKmbjxkGHDtC8ebyau8UW6bKsumrcLq5XL7j+eujRI57kJ0l5yH2SJalYTZwI++4bT8h77jnYYIPUieL4xS23xCznnBMXEvbvDyuumDqZJP2EV5IlqRhNnRqvIM+eDcOH50dBrurss6FfP3j2WejYMZZlScojlmRJKjZz5kCXLvDvf8OQIWlHLBbn+OPjwSOvvAKdO8N336VOJEn/z5IsScXmjDPgxRfhrrtgp51Sp1m8ww6DRx6Bl1+G3/wGZs5MnUiSAEuyJBWXu++Gvn3jOEOPHqnT1EzXrnDPPXFruq5dXcwnKS9YkiWpWLzzDpx6Kuy9N1x1Veo0S+eoo+C222DwYDjuOI+xlpScu1tIUjGYOTMe/7zSSvDgg7k9ZrqunHRSXHDYuzdsuCFcfnnqRJJKmCVZkorBuefCv/4VF+rV5lHTte2SS2DSJOjTB1q3hqOPTp1IUomyJEtSoSsvj3sP/+EP0KlT6jTLJwT4619hyhTo2RNatYL27VOnklSCQpaHc19lZWVZRUVF6hiSlP9mzICtt4ZGjeDNN2GFFVInyo1vv4VddoH//AcqKuL4hSTlWAhhbJZlZdU958I9SSpkF10EkyfHXS2KpSADNGsGTz0F8+fHreHcQ1lSHbMkS1KhGjUqjlmccgrsumvqNLn3i1/Ew0b+9S849lh3vJBUp5ZYkkMI94QQvgwhjF/E8+eEEN6svI0PIcwLIaxW+dykEMK/Kp9zfkKScmXu3LgbxPrrw5VXpk5Tezp2jH++xx+Ha69NnUZSCanJleT7gI6LejLLsmuzLGubZVlb4ALghSzLvqrykj0rn6923kOStAxuuy1eYb3xxrjtWzE791w49FC44AJ47rnUaSSViCWW5CzLXgS+WtLrKh0OPLpciSRJi/fFF3GrtA4d4OCDU6epfSHEmetNNoHu3eHLL1MnklQCcjaTHEJYkXjFeUCVhzNgeAhhbAihZ65+liSVtPPOgx9+iPPIIaROUzdWWgmeeAK+/hqOOCIu6JOkWpTLhXsHAKMWGrX4VZZl2wH7AaeEEHZb1BeHEHqGECpCCBVTp07NYSxJKiKvvw733w9nnAG//GXqNHVrm23g5pth+PDCO3ZbUsHJZUnuxkKjFlmWfVZ5/yXwd6Ddor44y7J+WZaVZVlW1rJlyxzGkqQict550KIFXHhh6iRpHH98PH77kkvgn/9MnUZSEctJSQ4hNAd2BwZWeaxpCGHlBR8DHYBqd8iQJNXA8OHxdL2LL4bmzVOnSWPBiXwbbQQ9esTxC0mqBTXZAu5R4BVg0xDClBDCsSGEE0MIJ1Z52cHA8CzLqu72vibwUgjhLeA1YEiWZc/mMrwklYz58+NV5Nat49ZvpaxZM3j4YfjsM+jVK3UaSUWqwZJekGXZ4TV4zX3EreKqPjYRaLOswSRJVTzySDx2+uGHoXHj1GnSa9cOLr0UeveGAw6Aww5LnUhSkQlZHp5gVFZWllVUePaIJAEwe3ZcpNeiBYwZA/U8LBWIB6rsuiu8/37cM3rddVMnklRgQghjF3WWh/+llaR8d999MHky/OlPFuSqGjSABx+EWbPg6KPdFk5STvlfW0nKZ7Nnw5//DDvuCPvumzpN/tlkE7jhBhgxAm69NXUaSUXEkixJ+eyBB+JV5EsvLZ2DQ5ZWz56w//5xYeOECanTSCoSlmRJyldz5sQRi7Iy2G+/1GnyVwjQrx+ssAIce6xjF5JywpIsSfnqoYdg0qS4g4NXkRdv7bXhpptg1CjHLiTlhLtbSFI+mjsXNtssHhpSUWFJroksg86d4YUX4m4XG22UOpGkPOfuFpJUaAYMgI8+iqfrWZBrZsHYRYMGcNxxsTRL0jKyJEtSvskyuPbauHPDgQemTlNY1lsPrrsOnnsuFmZJWkaWZEnKNy+8AGPHwllnQf36qdMUnuOOg/bt4Zxz4JNPUqeRVKAsyZKUb669Flq2hCOPTJ2kMIUAd94Zd7k44QTHLiQtE0uyJOWTt9+GZ56BU0+NW5pp2Wy4IVx5JTz7LDz+eOo0kgqQJVmS8sl118VyfNJJqZMUvpNPhh12gD/8AaZPT51GUoGxJEtSvvj8c3j4YTjmGFh99dRpCl/9+nHx3rRpcP75qdNIKjCWZEnKF3/9azxl7w9/SJ2keLRtC2ecEcvySy+lTiOpgFiSJSkfzJ4Nd9wBHTvGrd+UO5ddBq1axUV8s2enTiOpQFiSJSkf/P3vcdzi1FNTJyk+TZvCbbfBO+/EnUMkqQYsyZKUD265JR6j3LFj6iTFqXNn6NoVrrgCPvwwdRpJBcCSLEmpvfEGjBoFp5wC9fzPcq256SZo3BhOPNG9kyUtkf81lqTU+vaFFVeEo49OnaS4rbMOXHUVjBwJjzySOo2kPGdJlqSUvvoqbvv2u9/BqqumTlP8TjgB2rWLR35/803qNJLymCVZklJ66CGYOdPDQ+pKvXpxEd+XX8Kll6ZOIymPWZIlKZUsgzvvhLKyuJ+v6sb228e55FtvhbfeSp1GUp6yJEtSKq+9BuPHw/HHp05Sev70J1httXh09fz5qdNIykOWZElK5a674oK9bt1SJyk9q64K11wDL78MDzyQOo2kPGRJlqQU/vc/ePTRWJCbNUudpjQddRTssgucey5Mn546jaQ8Y0mWpBQefxy++w6OOy51ktK1YBHftGlw0UWp00jKM5ZkSUrhzjthyy1hp51SJyltbdpAr15wxx0wdmzqNJLyiCVZkurauHFx0d5xx0EIqdOoTx9YYw0X8Un6CUuyJNW1u+6CRo3giCNSJxFA8+Zw3XXxjctdd6VOIylPWJIlqS7Nnh1P2Dv4YGjRInUaLdCjB+y+O1xwAfz3v6nTSMoDlmRJqkvPPBOPoj7qqNRJVFUI0LdvPKr6ggtSp5GUByzJklSXHngA1lwT9tkndRItbMst4fTT4e674+iFpJJmSZakujJtGgweDN27Q4MGqdOoOr17w1prwSmnwLx5qdNISsiSLEl15YknYM4cOPLI1Em0KCuvDNdfDxUVLuKTSpwlWZLqygMPwNZbx715lb+6dYuL+C680EV8UgmzJEtSXXj/fRg9Om775t7I+S0EuPXWuIjvwgtTp5GUiCVZkurCQw/FY5B79EidRDWx1Vbwhz/EkYsxY1KnkZSAJVmSatv8+fDgg7D33rDOOqnTqKZ69447kZx8sov4pBJkSZak2vbSSzBpkifsFZpmzeJJfBUVcVs4SSXFkixJte2RR2DFFeMpeyos3bvDbrvFA0amTUudRlIdsiRLUm2aMwf694eDDoKmTVOn0dKquojvootSp5FUhyzJklSbRo6MVyC7dUudRMtq663h1FOhX784eiGpJFiSJak2PfYYNG8O++6bOomWx2WXwRprxJP45s9PnUZSHbAkS1JtmTkT/v53+O1voXHj1Gm0PJo3j4v4XnsN7rkndRpJdcCSLEm15dln4dtv4bDDUidRLvToAb/+NZx/Pnz1Veo0kmqZJVmSastjj8Hqq8Nee6VOolxYsIjv669dxCeVAEuyJNWG776Dp5+GLl2gYcPUaZQr22wDvXrBX/8KY8emTiOpFi2xJIcQ7gkhfBlCGL+I5/cIIXwTQniz8nZplec6hhAmhBA+DCGcn8vgkpTXnn4avv/eXS2K0eWXu4hPKgE1uZJ8H9BxCa/5Z5ZlbStvfQBCCPWBvsB+wBbA4SGELZYnrCQVjMcei0dQ77pr6iTKtebN4Zpr4NVX4d57U6eRVEuWWJKzLHsRWJYVCu2AD7Msm5hl2WzgMeCgZfg+klRYvvkGhg6FQw+F+vVTp1FtOOKI+AbIRXxS0crVTPLOIYS3QghDQwhbVj62LvBplddMqXxMkorb00/D7NnQtWvqJKotCxbxffUVXHxx6jSSakEuSvLrwAZZlrUBbgGeqnw8VPPabFHfJITQM4RQEUKomDp1ag5iSVIiAwbEUYuddkqdRLWpTZs4l3zHHfD666nTSMqx5S7JWZZ9m2XZjMqPnwEahhBWJ145Xr/KS9cDPlvM9+mXZVlZlmVlLVu2XN5YkpTGjBlxf+SDD4Z6biBU9Pr0gZYtXcQnFaHl/i94CGGtEEKo/Lhd5fecBowBNgkhbBhCaAR0AwYt78+TpLw2dGg8ae+QQ1InUV1YZZW4iG/0aLjvvtRpJOVQTbaAexR4Bdg0hDAlhHBsCOHEEMKJlS/pAowPIbwF/AXolkVzgV7AMOBd4Iksy96unT+GJOWJAQPiASK//nXqJKorRxwBu+wC550H06enTiMpR0KWLXJMOJmysrKsoqIidQxJWjozZ8ZfvXfrBnfemTqN6tKbb8L228NJJ8UFfZIKQghhbJZlZdU958CcJOXKiBFxJtlRi9LTti2cfDLcfju88UbqNJJywJIsSbkyYEA8aGKvvVInUQpXXAEtWriITyoSlmRJyoU5c2DQIDjwQGjUKHUapbBgEd8rr8ADD6ROI2k5WZIlKReeey4u2nLUorQdeSTsvDOce66L+KQCZ0mWpFwYMACaNoUOHVInUUr16kHfvjBtGlx6aeo0kpaDJVmSlte8efDUU9C5M6ywQuo0Sm3bbeMuF7fdFne9kFSQLMmStLxeegm+/NJRC/3IRXxSwbMkS9LyGjgwLtbbb7/USZQvVl0VrroKXn4ZHnwwdRpJy8CSLEnLI8virhbt28PKK6dOo3zy+9/DTjvFRXxff506jaSlZEmWpOXxzjvw0Udx6zepqgWL+KZOdRGfVIAsyZK0PAYNivcHHJA2h/LTdtvFRXx9+8LYsanTSFoKlmRJWh6DBkFZGay7buokyld/+hOssQb07Alz56ZOI6mGLMmStKw+/xxefdVRCy3eKqvAX/4Cr78e7yUVBEuyJC2rwYPjwr2DDkqdRPmuSxfYf3+45BKYNCl1Gkk1YEmWpGU1aBBssAFsvXXqJMp3IcS55BDg5JPjmytJec2SLEnL4vvvYcSIOGoRQuo0KgStWsX55KFD4YknUqeRtASWZElaFiNGwMyZziNr6fTqFRd6nnYaTJ+eOo2kxbAkS9KyGDQImjeH3XdPnUSFpH59uPNOmDYtHjIiKW9ZkiVpac2bB08/HY+hbtgwdRoVmrZt4ayz4K674IUXUqeRtAiWZElaWq+9Fk9Rc9RCy6p3b9hwQzjhhDi2IynvWJIlaWkNHAgNGsQrydKyWHFFuP12mDAB/vzn1GkkVcOSLElLa9CgOIu8yiqpk6iQ7bsv/O53cOWV8NZbqdNIWoglWZKWxocfwrvvOmqh3LjpJmjRAo4+GubMSZ1GUhWWZElaGkOGxPv990+bQ8WhRYs4dvHGG3D11anTSKrCkixJS2PIENhsM9hoo9RJVCwOPhgOOwz69IHx41OnkVTJkixJNTVjRtyyq3Pn1ElUbG65Jc64//73MHdu6jSSsCRLUs394x8wezZ06pQ6iYpNy5bQty+MHQvXXZc6jSQsyZJUc0OGwMorw667pk6iYtS1KxxySNxD+Z13UqeRSp4lWZJqIsvgmWdgn32gUaPUaVSs+vaNb8SOOSae7CgpGUuyJNXE+PEwZYqjFqpda64Z55NffRWuuSZ1GqmkWZIlqSaeeSbee8qealu3bnDooXDppXFrOElJWJIlqSaGDIFtt4V11kmdRMUuhLh3csuW8US+mTNTJ5JKkiVZkpZk+nR4+WVHLVR3VlsN7r03LuC78MLUaaSSZEmWpCUZMSIuorIkqy7tuy/06gU33ggjR6ZOI5UcS7IkLcmQIfHK3o47pk6iUnP11bDppvGQka+/Tp1GKimWZElanPnzYehQ6NgR6tdPnUalZsUV4aGH4PPP41VlSXXGkixJizN2LEyd6qiF0ikriztdPPwwPP546jRSybAkS9LiDBkSdxvo2DF1EpWyCy6AnXaCE06ASZNSp5FKgiVZkhbnmWdiOWnRInUSlbIGDeCRR+LJj4cfDnPmpE4kFT1LsiQtyhdfwJgxjlooP2y4Idx1F4weHccvJNUqS7IkLcqwYfG+c+e0OaQFunaFnj3hqqtg+PDUaaSiZkmWpEUZMgTWXhvatk2dRPrRjTfCllvCEUfEXS8k1QpLsiRVZ968eIjIvvvGhXtSvlhxRXjsMfj2WzjyyLhNoaScsyRLUnXGjInHUburhfLRVlvBzTfHN3JXX506jVSULMmSVJ1hw+IV5L33Tp1Eqt7xx8Nhh8HFF3tstVQLLMmSVJ1hw2CHHdz6TfkrhLjbxaabxm3hpkxJnUgqKpZkSVrY9Onw6qtxHlnKZyutBH/7G/zwQ9z5Yvbs1ImkomFJlqSFjRwZF0NZklUINtsM7r037p985pmp00hFw5IsSQsbNgyaN4cdd0ydRKqZLl1iQe7bFx5+OHUaqSgssSSHEO4JIXwZQhi/iOd7hBDGVd5eDiG0qfLcpBDCv0IIb4YQKnIZXJJqRZbFkty+fTwKWCoUV10Fv/51PGxk3LjUaaSCV5MryfcBi9sD6WNg9yzLtgGuAPot9PyeWZa1zbKsbNkiSlIdeu89+PRTRy1UeBo2hCeeiL8FOeggmDo1dSKpoC2xJGdZ9iLw1WKefznLsumVn44G1stRNkmqewuOorYkqxCttRY89VQ8ia9LFxfyScsh1zPJxwJDq3yeAcNDCGNDCD1z/LMkKfeGDYtbam2wQeok0rJp1w7uvhtefBFOOSWOEElaajkbuAsh7EksybtWefhXWZZ9FkJYAxgRQniv8sp0dV/fE+gJ0KpVq1zFkqSamzkTXnghHtIgFbLu3eHtt+HPf4att4bTTkudSCo4ObmSHELYBrgLOCjLsmkLHs+y7LPK+y+BvwPtFvU9sizrl2VZWZZlZS1btsxFLElaOv/8Z9xv1lELFYMrroizyWecAcOHp04jFZzlLskhhFbA34Ajsix7v8rjTUMIKy/4GOgAVLtDhiTlhWHDoFEj2H331Emk5VevHjz0EGy1FRx6KLz7bupEUkGpyRZwjwKvAJuGEKaEEI4NIZwYQjix8iWXAi2A2xba6m1N4KUQwlvAa8CQLMuerYU/gyTlxrBhcQutpk1TJ5FyY6WVYOBAaNwY9tsvLuiTVCNLnEnOsuzwJTx/HHBcNY9PBNr8/CskKQ/9+98wfjwceWTqJFJutW4NQ4bE35B07hzn7ldaKXUqKe954p4kwY8zm84jqxiVlcGTT8Jbb0HXrjBnTupEUt6zJEsSxFGLtdeOOwFIxahTJ7j9dnj2WTjpJLeGk5bAM1clad48GDECDjgAQkidRqo9xx8fT5S84gpYd124/PLUiaS8ZUmWpLFj4auvHLVQabj88jiD36cPrLJK3CJO0s9YkiVp2LB4BXmffVInkWpfCNCvH3z7LZx5JjRrBscemzqVlHcsyZI0bBhsvz2svnrqJFLdqF8fHn4YZsyAnj1jUe7aNXUqKa+4cE9Safv2Wxg9Gjp0SJ1EqluNGsGAAbDLLtCjBwwdmjqRlFcsyZJK2/PPx4V7jlqoFK24IgweHHd1Ofjg+FsVSYAlWVKpKy+PRWHnnVMnkdJo3jzuE7755nDQQXGLOEmWZEklbsQI2G23eGyvVKpatIhvGLfYIhblZ55JnUhKzpIsqXRNmQLvveeohQQ/FuWttoqjF0OGpE4kJWVJllS6ysvj/d57p80h5YvVVov/XiyYUe7fP3UiKRlLsqTSNWIErLGGR1FLVa26aizK7drBYYfBnXemTiQlYUmWVJqyLBaBvff2KGppYausEhfz7btv3Ef56qtTJ5LqnCVZUmn617/gyy+dR5YWZcUVYeBA6N4dzj8fzj03vrmUSoQn7kkqTc4jS0vWsCE8+GAcwbj2WvjsM7j7bneDUUmwJEsqTSNGwGabwXrrpU4i5bd69eCWW2DtteHii+GTT+Dvf4+7YUhFzHELSaVn1ix44QVHLaSaCgEuuggeewxeew122gnefz91KqlWWZIllZ5XXoEffnDUQlpahx0G//gHfP11PKXy+edTJ5JqjSVZUukZMQLq14c99kidRCo8u+wCo0fH7RP33htuvtkFfSpKlmRJpae8PP66uFmz1EmkwrTxxrEod+4Mp58ORxwB33+fOpWUU5ZkSaVl+nSoqHDUQlpezZvHBXxXXAGPPAK/+hV8/HHqVFLOWJIllZZ//APmz3fRnpQL9erFHS8GD4ZJk2C77WDAgNSppJywJEsqLeXlsPLK8chdSbnRqVP8Dc0vfgFdusCJJ8bFsVIBsyRLKi0jRsQFew0bpk4iFZeNN4ZRo+Ccc+Cvf4UddoDx41OnkpaZJVlS6fj4Y/joI+eRpdrSqBFccw0MGwZTp8aifMMNMG9e6mTSUrMkSyodC46idh5Zql0dOsC4cfHftbPOgt128/ARFRxLsqTSMWIErLtuPI5aUu1ac00YOBAeeADeeQfatIEbb/SqsgqGJVlSaZg/H0aOjKMWIaROI5WGEOIeym+/Hf/dO/PMeFLf2LGpk0lLZEmWVBreeAO++spRCymFddaBQYPg4Yfhk0/i7jKnnhqPt5bylCVZUmlYMI/cvn3aHFKpCgG6d4f33oNTToHbboujTw88EH/TI+UZS7Kk0jBiBGy9Nay1VuokUmlbZRX4y19gzBjYYAM46qi4C8Zzz6VOJv2EJVlS8fvhB3jpJUctpHyy3XbwyitxBOO//4W99oIDD4R3302dTAIsyZJKwUsvwaxZ7o8s5Zt69eIIxoQJcNVV8MILsNVWcOSRbhmn5CzJkorfiBHxhL3ddkudRFJ1mjSB886DDz+M+yoPGACbb25ZVlKWZEnFr7wcdtkFmjZNnUTS4rRsGU/s+/jjuF1c//6xLPfoEXeokeqQJVlScZs6Nf7l6jyyVDjWWAOuvfbHsjxoUJxh3nNPGDzY3TBUJyzJkorbyJHx3pIsFZ4114xlecqUeP/hh3DAAfHqct++7rOsWmVJllTcysvjllPbb586iaRl1bw5nH02TJwIjzwCzZpBr17xkJLf/x5GjYIsS51SRcaSLKl4ZVlctLfXXlC/fuo0kpZXw4Zw+OFxj+WKiriw729/g113jbtiXHstfPpp6pQqEpZkScXrww/jEbiesicVn+23hzvugM8+g7vuileXzz0XWrWKO9nccUfcf1laRpZkScXLeWSp+K20Ehx7bDyY5IMP4IorYjk+6SRYe+24P/ott8DkyamTqsCELA9neMrKyrKKiorUMSQVui5d4q9lJ02CEFKnkVRXsgzGjYPHHoOBA388xa9Nm3iq3/77xyvRjmGVvBDC2CzLyqp9zpIsqSjNmxf3XD34YLj77tRpJKX0wQdxG7lBg+IJnPPnxwW9e+wRrzS3bw+bbuqb6RK0uJLcoK7DSFKdeOMNmD7do6glwSabxJP8zjorjmKUl8dxrPJyeOqp+Jp11oFf/xp23jne2raFRo3S5lZSlmRJxam8PN7vtVfaHJLyy+qrQ7du8QZxW7mRI+Pt5Zfh8cfj402aQFlZLMxlZbE0b7yxIxolxHELScVp773hyy/jXKIk1dS//x0XAb78crx//XWYPTs+t+KKsM02sTC3bQtbbgmbbRaLtwqSM8mSSssPP8Cqq8LJJ8MNN6ROI6mQzZoVF/69+eZPb9988+NrVlstzjRvtlm833RT2Ggj2GCDeBCK8pYzyZJKy8svx7/YnEeWtLwaN/7xyvECWRa3lHvnHZgwId7eew+GDoV77/3p16+yCrRuHQtz69bxtu66cXu6tdaK902b1uEfSDVVo5IcQrgH2B/4Msuyrap5PgA3A52A74HfZ1n2euVzRwEXV770j1mW3Z+L4JK0SOXl0KBBPFBAknIthB8Lb6dOP33um2/g/ffj1pMLbpMnx8ONysvhu2Gs9E0AACAASURBVO9+/v1WWunH0rzg1qJFvEK98K1Fi3h12tnoWlejcYsQwm7ADOCBRZTkTsCpxJK8I3BzlmU7hhBWAyqAMiADxgLbZ1k2fXE/z3ELSctlhx1ghRXgxRdTJ5GkH2UZfPVVnHv+/HP4z39+er/g4y+++Ok4x8JCiEW5WbNYsFdeOd4W9/EKK8TFiE2a/PTjhW8rrBCvnpfIdnjLPW6RZdmLIYTWi3nJQcQCnQGjQwirhBDWBvYARmRZ9lVlkBFAR+DRmsevG88++yyff/556hiSllOjGTM4fOxY3jzwQN66777UcSRp8dZcM94WEubNo9F339G4yq3Rd9/ReMaM//+84cyZNJg5k4YzZtDwv/+l4cyZ8bFZs2g4cyb1585d5ljzGjRgXoMGZPXrM79+febXq/f/H2f16sXHqn7eoMFPXlP1dVkIEAJZvXpki7j/fo89KLvoouX5J5lzuZpJXhf4tMrnUyofW9TjPxNC6An0BGjVqlWOYkkqNWu/9x4hy/jPllumjiJJyyyrX59ZzZoxq1mzZf4e9ebOjSV65kzqz5mzyFuD6h6fPZt68+ZRb948wvz51Js7l3rz5xMqH/v/x6t8XH/uXOrNmhVfU/lcmD8/3rLs/z8my6hXeb/guY832CCH//RyI1clubpr8tliHv/5g1nWD+gHcdwiR7lqrGPHjnX9IyXVhldfhZVWotNll0HDhqnTSJJqYLPUAapRL0ffZwqwfpXP1wM+W8zjklQ7Ro6MR81akCVJyyFXJXkQcGSIdgK+ybLsP8AwoEMIYdUQwqpAh8rHJCn3Jk+GDz5w6zdJ0nKr6RZwjxIX4a0eQpgC9AYaAmRZdgfwDHFniw+JW8AdXfncVyGEK4Axld+qz4JFfJKUcyNHxntLsiRpOdV0d4vDl/B8BpyyiOfuAe5Z+miStJTKy+P+oltskTqJJKnA5WrcQpLSmj8/luT27Utmf09JUu2xJEsqDuPHw9SpjlpIknLCkiypOCyYR27fPm0OSVJRsCRLKg7l5bDpprD++kt+rSRJS2BJllT4Zs+GF15w1EKSlDOWZEmF79VX4bvvLMmSpJyxJEsqfOXlUK9ePGlPkqQcsCRLKnzl5VBWBquskjqJJKlIWJIlFbZvv43jFo5aSJJyyJIsqbC9+CLMm2dJliTllCV5gQcegEceiX/ZSioc5eWwwgqw886pk0iSiogleYH774cePaBNG/jb3yDLUieSVBPl5fDrX0OTJqmTSJKKiCV5gREj4LHHYO5cOOQQ2H57GDLEsizls//8B95+21P2JEk5Z0leoF49OOwwGD8+XlX++mvYf3/YbTeoqEidTlJ1FhxF7TyyJCnHLMkLa9AAjjwSJkyAO+6A99+HHXaA3/8ePvssdTpJVY0cCautBm3bpk4iSSoyluRFadgQTjgBPvgAzjsPHn0UNtkE/vQnmDUrdTpJWRbnkdu3j78JkiQph/ybZUmaNYOrroJ334WOHeHii+NVqxdfTJ1MKm3vvw9TpjhqIUmqFZbkmtpoIxgwAJ55BmbOhN13h2OPhWnTUieTSlN5eby3JEuSaoEleWntt19cTX/eeXFv5c02g/79U6eSSk95ObRuHd/ASpKUY5bkZbHiinEE4/XX41/SXbvGPZanT0+dTCoNc+fCc895FVmSVGssyctj663h5ZfhssvgiSdgq61g2LDUqaTiN3YsfPONJVmSVGssycurYUPo3RtGj4ZVVomL+047zR0wpNq0YH/kvfZKm0OSVLQsybmy/fbx6tYf/gC33AK/+hV89FHqVFJxKi+Pu8y0bJk6iSSpSFmSc6lJE7jpJnjqqViQt9sOnnwydSqpuHz/PYwa5aiFJKlWWZJrw0EHwZtvwhZbwKGHQq9eMGdO6lRScXjpJZg9Ox4iIklSLbEk15YNNogHjpx1FvTtG696ffll6lRS4Ssvj2sBfv3r1EkkSUXMklybGjaE666Dhx+G116DsrK4bZykZTd8eJz5b9o0dRJJUhGzJNeF7t3jDCXEv9wfeSRtHqlQffEFvPUWdOiQOokkqchZkuvKdttBRQW0axcPHrniCsiy1KmkwrLgKGpLsiSpllmS69Iaa8CIEXDkkXDppXDccS7ok5bGiBHQogVsu23qJJKkItcgdYCS06gR3HdfPM66Tx+YMiVuE9esWepkUn7LsjiPvPfeUM/395Kk2uXfNCmEAJdfDnffHU8O2203+M9/UqeS8tvbb8d/T/bZJ3USSVIJsCSndMwxMGQIfPhh3M5q0qTUiaT8NWJEvLckS5LqgCU5tX33jYuRpk2LRfm991InkvLT8OGw6abQqlXqJJKkEmBJzgc77QQvvBAX8e22G7zxRupEUn6ZNSv+O+KuFpKkOmJJzhfbbAP//CessALsuSe88krqRFL+GDUKfvjBkixJqjOW5HyyySaxKK+xRhzDGD06dSIpP4wYAQ0awO67p04iSSoRluR806oVPPfcj0X5tddSJ5LSGz4cdtkFVl45dRJJUomwJOejddeNRXn11eOvlysqUieS0pk6FV5/3V0tJEl1ypKcr9ZfPxblVVeN5WDs2NSJpDRGjoz3ziNLkuqQJTmfLRi9aN48jl68+27qRFLdGz48vlncfvvUSSRJJcSSnO9at477KDdoEK8oT56cOpFUd7IsLtpr3x7q10+dRpJUQizJheAXv4hX02bMiL9y/vLL1ImkuvHeezBliqMWkqQ6Z0kuFNtsA4MHw6efwn77wbffpk4k1b7hw+O9i/YkSXXMklxIdt0V+veHcePgwANh5szUiaTaNWJE3D+8devUSSRJJcaSXGg6dYL7749H9B59NMyfnzqRVDtmz4bnn3fUQpKURIPUAbQMuneHTz6BCy6AjTaCP/0pdSIp9155Bb77zlELSVISluRCdd55MHEi/PnPsOGGcNxxqRNJuTV8eNzRYs89UyeRJJWgGo1bhBA6hhAmhBA+DCGcX83zN4YQ3qy8vR9C+LrKc/OqPDcol+FLWgjQt2/cP/nEE2HYsNSJpNwaPhx22gmaNUudRJJUgpZYkkMI9YG+wH7AFsDhIYQtqr4my7Izsixrm2VZW+AW4G9Vnv5hwXNZlh2Yw+xq2BCeeAK23BK6dIF//St1Iik3pk2Lp0w6jyxJSqQmV5LbAR9mWTYxy7LZwGPAQYt5/eHAo7kIpxpo1gyGDIGVV4aDDoL//jd1Imn5lZfHg0ScR5YkJVKTkrwu8GmVz6dUPvYzIYQNgA2Bf1R5uEkIoSKEMDqE8JtlTqpFW289+Pvf4d//hkMPhTlzUieSls+zz8ajqNu1S51EklSialKSQzWPZYt4bTegf5Zl86o81irLsjKgO3BTCGHjan9ICD0ry3TF1KlTaxBLP7HjjtCvHzz3HJx5Zuo00rLLsliSO3TwKGpJUjI1KclTgPWrfL4e8NkiXtuNhUYtsiz7rPJ+IvA8sG11X5hlWb8sy8qyLCtr2bJlDWLpZ446Cs44A269Fe66K3UaadmMGweffw4dO6ZOIkkqYTUpyWOATUIIG4YQGhGL8M92qQghbAqsCrxS5bFVQwiNKz9eHfgV8E4ugmsRrrkmznGefDKMGpU6jbT0nn023u+7b9ockqSStsSSnGXZXKAXMAx4F3giy7K3Qwh9QghVd6s4HHgsy7KqoxibAxUhhLeA54CrsiyzJNemBg3gscegVSv47W/jnLJUSJ59Ftq0gbXXTp1EklTCwk87bX4oKyvLKioqUscobO+8Exc9tW0b55QbNkydSFqyb7+FFi3g7LPhyitTp5EkFbkQwtjKtXM/U6PDRFSAttgiziWPGgXn/+z8Fyk//eMfMHeu88iSpOQsycWsWzc45RS44Qb429+W/HoptWefjXt+77xz6iSSpBJnSS52118fxy6OPho++CB1GmnRFmz91r49NGqUOo0kqcRZkotd48bx6OoGDeLR1d9/nzqRVL0JE2DyZEctJEl5wZJcCjbYAB56KO4/e8opqdNI1XPrN0lSHrEkl4r99oOLL4b77ouFWco3Q4fCZptB69apk0iSZEkuKb17w667wkknwUcfpU4j/ej77+GFF+KbOUmS8oAluZQ0aAAPPxzvDz8cZs9OnUiKXngBZs1yHlmSlDcsyaWmVSu4804YMwYuvTR1Gil69llYYQXYbbfUSSRJAizJpalLF+jZE66+GkaMSJ1GiiV5jz2gSZPUSSRJAizJpevGG2HzzeHII2Hq1NRpVMomToT333fUQpKUVyzJpWrFFeHRR2H69HjQSJalTqRSNWxYvLckS5LyiCW5lLVpA9dcA0OGxDllKYWhQ2HDDWGTTVInkSTp/1mSS12vXvEY4DPPdFs41b0ffoCRI6FTJwghdRpJkv6fJbnU1asH994bt4U76iiYNy91IpWS55+PeyTvv3/qJJIk/YQlWbD++nDrrTBqFFx/feo0KiVDhsT5+D32SJ1EkqSfsCQr6tEDDjkELrkExo1LnUalIMtiSd57b7d+kyTlHUuyohDg9tth1VXhiCPi6WdSbXrnHZg0CTp3Tp1EkqSfsSTrRy1bxl0uxo2Dyy9PnUbFbsiQeN+pU9ockiRVw5KsnzrgADj22Hga3yuvpE6jYjZ4MLRtC+utlzqJJEk/Y0nWz91wQywuRx8NM2emTqNi9NVX8PLLjlpIkvKWJVk/16xZHLuYMAH69EmdRsVo2LC43aBbv0mS8pQlWdXr0AGOOSaeyDd2bOo0KjZDhsDqq8MOO6ROIklStSzJWrTrr4c11ohlefbs1GlULObNi0dRd+oE9eunTiNJUrUsyVq0VVaBO+6Iu11cdVXqNCoWo0fHmWTnkSVJecySrMU78EDo3h3++EcYPz51GhWDIUPiMegdOqROIknSIlmStWQ33xyvKh9zDMydmzqNCt3gwbDrrvH/U5Ik5SlLspZs9dXh1lthzJi4PZy0rD75BP71L3e1kCTlPUuyaqZrVzj4YOjdGz76KHUaFapnnon3ziNLkvKcJVk1E0K8mtyoEZx4ImRZ6kQqRIMGwcYbw6abpk4iSdJiWZJVc+usA1deCeXl8PDDqdOo0PzvfzByJBx0UHzTJUlSHrMka+mceCLstBOccQZMm5Y6jQrJsGFxv+2DDkqdRJKkJbIka+nUqwf9+sHXX8M556ROo0IycCC0aAG77JI6iSRJS2RJ1tLbems4+2y49154/vnUaVQI5syJW78dcEDcI1mSpDxnSdayufRS2GgjOOEEmDkzdRrlu3/+M/72wVELSVKBsCRr2aywQjyy+v3342I+aXEGDoQmTWCffVInkSSpRizJWnb77AM9esSS/O67qdMoX2VZLMn77ANNm6ZOI0lSjViStXxuuAFWXjmOXbh3sqrz1lsweTL85jepk0iSVGOWZC2fNdaAa66JM6cPPJA6jfLRwIFxX2SPopYkFRBLspbf0UfDzjvHLeGmT0+dRvlm4MC47dsaa6ROIklSjVmStfzq1YPbb4+Hi1x0Ueo0yieffAJvvOGohSSp4FiSlRtt2sCpp8YdL8aMSZ1G+WLgwHjv1m+SpAJjSVbu9OkDa60FJ58M8+alTqN8MHAgbL45bLJJ6iSSJC0VS7Jyp1kzuP56qKiIR1ertE2bFk9kdNRCklSALMnKrW7dYM894cIL4csvU6dRSgMHxt8odOmSOokkSUvNkqzcCgH69oXvvoNzz02dRikNGAAbbgjbbps6iSRJS82SrNzbfHM46yy4//64f7JKz9dfw4gR8SpyCKnTSJK01CzJqh0XXwytWsVFfHPmpE6juvb00/F/d0ctJEkFypKs2tG0Kdx8M4wfD3/5S+o0qmv9+8P668MOO6ROIknSMqlRSQ4hdAwhTAghfBhCOL+a538fQpgaQniz8nZcleeOCiF8UHk7KpfhlecOOgg6d4bLLoMpU1KnUV353/9g2DA45BBHLSRJBWuJJTmEUB/oC+wHbAEcHkLYopqXPp5lWdvK212VX7sa0BvYEWgH9A4hrJqz9MpvIcSryHPmwNlnp06jujJkCMya5aiFJKmg1eRKcjvgwyzLJmZZNht4DKjp8Vn7AiOyLPsqy7LpwAig47JFVUHaaCO44AJ4/HEYOTJ1GtWF/v1h7bVh551TJ5EkaZnVpCSvC3xa5fMplY8t7JAQwrgQQv8QwvpL+bUqZueeG8tyr14we3bqNKpN330HzzwTRy3queRBklS4avK3WHVDhdlCnz8NtM6ybBugHLh/Kb42vjCEniGEihBCxdSpU2sQSwVjhRXi2MV778FNN6VOo9o0dCj88IOjFpKkgleTkjwFWL/K5+sBn1V9QZZl07Ism1X56Z3A9jX92irfo1+WZWVZlpW1bNmyJtlVSDp3hgMPhD594NNPl/x6Fab+/WGNNWDXXVMnkSRpudSkJI8BNgkhbBhCaAR0AwZVfUEIYe0qnx4IvFv58TCgQwhh1coFex0qH1MpuummeEzxWWelTqLa8MMPcdHewQdD/fqp00iStFyWWJKzLJsL9CKW23eBJ7IsezuE0CeEcGDly04LIbwdQngLOA34feXXfgVcQSzaY4A+lY+pFG24IVx4ITz5ZDyNTcXlmWdgxgzo2jV1EkmSllvIsmpHhJMqKyvLKioqUsdQbZg5E7baKl5pHDcOGjdOnUi50qULjBoV98T2SrIkqQCEEMZmWVZW3XMuP1fdatIEbrkF3n8fbrwxdRrlyrffwuDBcOihFmRJUlGwJKvu7bcf/OY3cMUV8MknqdMoF556Kh4gcvjhqZNIkpQTlmSlceONkGVw5pmpkygXHn00zpzvuGPqJJIk5YQlWWm0bg0XXQQDBsAwNzwpaFOnxoWY3brFo8glSSoClmSlc/bZ8ItfwKmnxl/VqzA9+WTc2s9RC0lSEbEkK53GjeMivg8+gOuvT51Gy+rRR2HLLWHrrVMnkSQpZyzJSqtjR/jtb+GPf4TJk1On0dL65BN46SWvIkuSio4lWekt2ArujDPS5tDSe/zxeN+tW9ockiTlmCVZ6bVqBZdcAn//OwwdmjqNlsYjj0C7drDxxqmTSJKUU5Zk5Yczz4Rf/jIu4ps5M3Ua1cS4cfDmm3DEEamTSJKUc5Zk5YcFi/g++giuuy51GtXE/fdDw4bOI0uSipIlWfmjQwfo0gX+9CeYNCl1Gi3O3Lnw8MOw//7QokXqNJIk5ZwlWfnlhhugXj04/fTUSbQ4w4bBF1/AkUemTiJJUq2wJCu/rL8+XHopDBwIQ4akTqNFuf/+eAW5U6fUSSRJqhWWZOWfM86ATTeF005zEV8+mj49vonp3h0aNUqdRpKkWmFJVv5p1AhuvRUmToRrrkmdRgt7/HGYPRuOOip1EkmSao0lWflp773h0EPhyivh449Tp1FV998fj6HebrvUSSRJqjWWZOWv66+H+vVdxJdPJkyA0aPjVeQQUqeRJKnWWJKVv9ZbD3r3hkGDYPDg1GkEcN99cfeRHj1SJ5EkqVZZkpXf/vAH2HzzuIjvhx9Spyltc+bAvfdC586wzjqp00iSVKssycpvCxbxffyxi/hSGzw47o3cs2fqJJIk1TpLsvLfXntBt25xEd/EianTlK5+/WDddaFjx9RJJEmqdZZkFYbrroOGDeP4here5MnxlL1jj4UGDVKnkSSp1lmSVRjWXRcuuyz+yv/pp1OnKT133x3vjzkmbQ5JkuqIJVmF47TTYIstXMRX1+bOhXvuiWMWG2yQOo0kSXXCkqzC0bAh9O0LkybBVVelTlM6hg6Ff/8bjj8+dRJJkuqMJVmFZY894PDD4eqr4cMPU6cpDbfdBmuvDfvvnzqJJEl1xpKswrNgEd9pp0GWpU5T3D74AJ59Fk48Mf4zlySpRFiSVXjWWQcuvzyOAQwalDpNcevbN5Zj90aWJJUYS7IK06mnwpZbxi3hvv8+dZriNGNGPGGva1dYa63UaSRJqlOWZBWmBYv4Jk+Oh4wo9x58EL79Nr4hkSSpxFiSVbh23x169IjHVX/wQeo0xSXL4nHg228PO+6YOo0kSXXOkqzCdu210LhxvNrpIr7cee45eOed+M81hNRpJEmqc5ZkFba114Y+feKRyU8+mTpN8bjxRlh9dTjssNRJJElKwpKswterVxwLOPVUmDYtdZrC98478fjvXr2gSZPUaSRJSsKSrMLXoAHcfTd89RWccUbqNIXv+uthhRXglFNSJ5EkKRlLsopDmzZw/vlxR4ahQ1OnKVyffRb/GR59dBy3kCSpRFmSVTwuvhi22AJOOCFuXaal95e/wLx5cOaZqZNIkpSUJVnFo3HjOHYxZUq8qqyl87//wR13wG9/CxtvnDqNJElJWZJVXHbaCU4/HW6/HV58MXWawtKvH3zzDZxzTuokkiQlF7I83Fu2rKwsq6ioSB1Dheq772CbbaB+fXjrrbgITYv3ww+w0Uaw2WZxj2RJkkpACGFslmVl1T3nlWQVn6ZN4c474yl8l16aOk1h6NcPPv8cevdOnUSSpLxgSVZx2msv6Nkzbmf20kup0+S3mTPh6qvjMd977JE6jSRJecGSrOJ1/fWw4YZw5JFxUZqqd+ed8J//eBVZkqQqLMkqXiutBPffD5MmwVlnpU6Tn2bOhKuugt128yqyJElVWJJV3HbdFc49N14tHTIkdZr8c/fd8QCR3r0hhNRpJEnKG+5uoeI3axa0awdffAHjx3uS3AIzZsAvfgGbbgrPP29JliSVHHe3UGlr3DgetTx9Opx4IuThG8Mkrr8+vnG45hoLsiRJC6lRSQ4hdAwhTAghfBhC+NlRZiGEM0MI74QQxoUQRoYQNqjy3LwQwpuVt0G5DC/V2DbbwB//CAMGwF13pU6T3uefw7XXQpcusOOOqdNIkpR3lliSQwj1gb7AfsAWwOEhhC0WetkbQFmWZdsA/YFrqjz3Q5ZlbStvB+Yot7T0zjoLOnSA006LYxel7PLL4xjKn/+cOokkSXmpJleS2wEfZlk2Mcuy2cBjwEFVX5Bl2XNZln1f+eloYL3cxpRyoF49eOABWGUVOPTQeDJfKZowIS5kPOEE2GST1GkkScpLNSnJ6wKfVvl8SuVji3IsMLTK501CCBUhhNEhhN8sQ0Ypd9ZcEx56CN57L15RLkXnnBOP6vY0QkmSFqkmJbm6FT3VrnwKIfwOKAOurfJwq8pVg92Bm0IIGy/ia3tWlumKqVOn1iCWtIzat4cLL4R77oFHHkmdpm4NGQJPPw0XXwxrrJE6jSRJeasmJXkKsH6Vz9cDPlv4RSGEvYGLgAOzLJu14PEsyz6rvJ8IPA9sW90PybKsX5ZlZVmWlbVs2bLGfwBpmVx2WdxD+YQT4vhBKZg5M14932wzOOOM1GkkScprNSnJY4BNQggbhhAaAd2An+xSEULYFvgrsSB/WeXxVUMIjSs/Xh34FfBOrsJLy6xBA3j0UWjSBH7zG/j229SJat+118LEiXDLLdCoUeo0kiTltSWW5CzL5gK9gGHAu8ATWZa9HULoE0JYsFvFtcBKwJMLbfW2OVARQngLeA64KssyS7Lyw3rrwZNPwgcfwFFHwfz5qRPVng8/jDtZdO0Ke++dOo0kSXnPE/ekm2+G00+HPn3gkktSp8m9+fNhr73gjTfg7bfjmwNJkrTYE/ca1HUYKe+cdhq8/jr07g3bbgv77586UW799a/wwgvxEBULsiRJNeKx1FIIcMcdsSB37w7jxqVOlDuTJ8O558YRi2OOSZ1GkqSCYUmWIO4bPHAgrLwydO4M//536kTLb968OGudZfHwkFDdbo6SJKk6lmRpgfXWi/sIf/11LMr/+1/qRMvnyivjmEXfvtC6deo0kiQVFEuyVFXbttC/P4wfH3eCmD07daJlM2pU3Au6e3c48sjUaSRJKjiWZGlh++4bF7sNGwa/+x3MnZs60dKZOjWW41at4PbbHbOQJGkZuLuFVJ1jj4VvvoGzzoIVV4xHWNcrgPeUc+bEK+BffAEvvQTNmqVOJElSQbIkS4ty5pkwY0bcGq5pU7j11vy/Knv66XEO+aGHoKzabR8lSVINWJKlxbnkEvjuO7jmGpg1K45h1K+fOlX1brkFbrsNzjkHevRInUaSpIJmSZYWJwS46ipo3BiuuCLuePHgg9CoUepkP/Xww/FQlIMOirtaSJKk5WJJlpYkhHhkdfPmcPbZsSg//njcUzkfDB4c90Pec0947LH8vdItSVIBKYCVSFKeOOuseCjH8OHwq1/Bxx+nTgR/+xscckjcuu6pp6BJk9SJJEkqCpZkaWkcdxwMHQqffgrt2sGLL6bLcu+9cSeL7baLxd2dLCRJyhlLsrS09tkHXn0VWrSA9u3jDPC8eXX38+fNiwsKjzkm/vzyclhttbr7+ZIklQBLsrQsfvlLGD0aDj4YLrwwltXJk2v/5371Fey/P/zxj3D00fD003F7OkmSlFOWZGlZrbJKXMB3770wdixssUXcCaO2jrIeOhTatIGRI+GOO+Duu+OuG5IkKecsydLyCAF+/3t4++14nPUFF8DWW8ct2XI1gjFxIhx2GHTqFOeOR42CE07I/4NNJEkqYJZkKRdatYo7TQwZEq/u/u538cryrbfGEYll8fbbcPzxsPnmcZu33r3h9ddhhx1ym12SJP2MJVnKpU6d4M03oX9/WGklOPVUWHvtOLt8220wYQLMn1/9186ZExcEXn11PFJ6q63iFeljj4UPPoDLLnO8QpKkOhKyLEud4WfKysqyioqK1DGk5ffGG3FmeeBA+OST+FiTJrDRRv/X3h29SFXGYRz/PoRRWpGghLiSodFNFxWLN2IXUWEl1WVBgRB0U2B0EXXZPxBdR3URRRJYEAWVkBlCpa5pZVuhYbQYWoSUIIb16+IMsZ0kiGb2HXe+Hxj2nJ2FeeDHMA/vvOcsLF/eHZ85AydOdBf+nTvX/c30dLfFYutWWLGiWXxJkhazJDNVHtlKCAAAA8RJREFUNX3e5yzJ0gKogqNHYfdumJ3t9hmfOgVnz8LSpd3t5Nat6y7M27SpW32WJEkj9W8l2X9LLS2EBNav7x6SJGnsuSdZkiRJ6rEkS5IkST2WZEmSJKnHkixJkiT1WJIlSZKkHkuyJEmS1GNJliRJknosyZIkSVKPJVmSJEnqsSRLkiRJPZZkSZIkqceSLEmSJPVYkiVJkqQeS7IkSZLUY0mWJEmSeizJkiRJUo8lWZIkSeqxJEuSJEk9qarWGf4hyY/Adw1eegXwU4PX1cJyzpPBOS9+zngyOOfJ0GrOV1fVyvM9MZYluZUk+6tqunUOjZZzngzOefFzxpPBOU+GcZyz2y0kSZKkHkuyJEmS1GNJ/rvnWgfQgnDOk8E5L37OeDI458kwdnN2T7IkSZLU40qyJEmS1GNJHkiyOcnXSY4kebJ1Hg1fkheTnEzyRessGo0ka5LsSjKb5HCSba0zafiSXJJkb5JDgzk/3TqTRiPJRUk+TfJW6ywajSTHknye5GCS/a3zzOd2C7o3IfANcBswB+wD7q+qL5sG01AluRk4DbxUVde3zqPhS7IKWFVVB5JcDswA9/peXlySBFhWVaeTLAH2ANuq6uPG0TRkSR4HpoErqmpL6zwaviTHgOmqGrt7YbuS3NkAHKmqb6vqN2A7cE/jTBqyqvoQ+Ll1Do1OVf1QVQcGx78Cs8Dqtqk0bNU5PThdMni44rPIJJkC7gKeb51Fk8mS3FkNfD/vfA4/WKULWpK1wI3AJ22TaBQGX8MfBE4CO6vKOS8+zwJPAH+0DqKRKuC9JDNJHm4dZj5Lcifn+Z2rEtIFKsllwA7gsar6pXUeDV9V/V5VNwBTwIYkbqFaRJJsAU5W1UzrLBq5jVV1E3AH8Mhga+RYsCR35oA1886ngOONskj6HwZ7VHcAr1TV663zaLSq6hTwAbC5cRQN10bg7sF+1e3ALUlebhtJo1BVxwc/TwJv0G2BHQuW5M4+4Nok1yS5GLgPeLNxJkn/0eCCrheA2ap6pnUejUaSlUmuHBxfCtwKfNU2lYapqp6qqqmqWkv3mfx+VT3QOJaGLMmywUXWJFkG3A6MzR2oLMlAVZ0DHgXepbvQ57WqOtw2lYYtyavAR8B1SeaSPNQ6k4ZuI/Ag3arTwcHjztahNHSrgF1JPqNb5NhZVd4iTLrwXAXsSXII2Au8XVXvNM70F28BJ0mSJPW4kixJkiT1WJIlSZKkHkuyJEmS1GNJliRJknosyZIkSVKPJVmSJEnqsSRLkiRJPZZkSZIkqedPx40rJjQ5UBcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "a,b = 0,5\n",
    "xd = np.linspace(a,b,1000)\n",
    "plt.plot([a,b],[1,1],c='grey')  # 1 line\n",
    "plt.plot(xd,dF(xd),c='red')  # plot the derivative\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Limitations of successive approximations method\n",
    "\n",
    "- Only stable solutions can be computed with SA  \n",
    "- Starting values matter for convergence!  \n",
    "- Also, only linear (slow) convergence to the solution when it does converge  \n",
    "\n",
    "\n",
    "But very easy to implement!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Multivariate example: platform market equilibrium\n",
    "\n",
    "- Platform is product or service with *network effects*, i.e. value (for firm and users) is generated from the interaction of users, and therefore depends on the size of the user base.  \n",
    "- Examples: game consoles, interface standards, social media, two-sided market platforms (Amazon, Ebay), peer-to-peer markets (Uber, AirBnB)  \n",
    "- Because of better experience of the users on larger platforms, platform markets tend to monopolization (*tipping markets*)  \n",
    "- Yet, there can be opposite forces as well: heterogeneous preferences or *congestion* in too large platforms  \n",
    "- Many different aspects (see further learning resources), we’ll look at a simple example, focusing on consumer choice and equilibrium  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### The model\n",
    "\n",
    "- $ m $ products, $ i=1,\\dots,m $  \n",
    "- unit mass of consumers of $ n $ types with different preferences over the product  \n",
    "- $ p_j $ are fractions of consumer types in the population, $ j=1,\\dots,n $  \n",
    "- utility function of consumers of type $ j \\in 1,\\dots,n $ from product $ i \\in 1,\\dots,m $ is given by  \n",
    "\n",
    "\n",
    "$$\n",
    "u_{ij} = c_{ij} + s_i,\n",
    "$$\n",
    "\n",
    "- $ c_{ij} $ are *valuations* of each product by each of the consumer type  \n",
    "- $ s_i $ are market shares of each product, increase utility when more people are using the same platform  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Choice probabilities\n",
    "\n",
    "Standard random utility framework with logit choice probabilities\n",
    "\n",
    "- value of choice $ U_{ij} = u_{ij} + \\epsilon[i] $ is effected by the random factors  \n",
    "- one $ \\epsilon[i] $ for each product choice, forming the vector $ \\epsilon = (\\epsilon[1],\\dots,\\epsilon[m]) $  \n",
    "- elements of $ \\epsilon $ are extreme value type I (EV1) random components of the utility, iid between consumers and products  \n",
    "\n",
    "\n",
    "$$\n",
    "P_{ij} = \\frac{\\exp(u_{ij})}{\\sum_{k=1}^m \\exp(u_{kj})}\n",
    "= \\frac{\\exp(u_{ij}-\\alpha)}{\\sum_{k=1}^m \\exp(u_{kj}-\\alpha)}, \\; \\forall \\alpha\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Market shares\n",
    "\n",
    "- assume choice is made once, and no changes are made afterwards  \n",
    "- under the assumption of unit mass of consumers we have  \n",
    "\n",
    "\n",
    "$$\n",
    "s_i = \\sum_{j=1}^n p_j P_{ij} = P \\cdot p\n",
    "$$\n",
    "\n",
    "- simple matrix product, assuming $ p $ is a column-vector  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Fixed point equation\n",
    "\n",
    "Combining the last three expressions:\n",
    "\n",
    "$$\n",
    "u_{ij} = c_{ij} + \\sum_{t=1}^n p_t \\frac{\\exp(u_{it})}{\\sum_{k=1}^m \\exp(u_{kt})}\n",
    "$$\n",
    "\n",
    "- $ mn $ fixed point equations  \n",
    "- fixed point in the ($ mn $)-dimensional space of $ u_{ij} $  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "class model:\n",
    "    '''Simple platform equilibrium model'''\n",
    "\n",
    "    def __init__(self,m=2,n=2):\n",
    "        '''Define default model parameters'''\n",
    "        self.m,self.n = m,n  # number of products and consumer types\n",
    "        self.c = np.random.uniform(size=(m,n))    # utilities (random uniform)\n",
    "        self.p = np.random.dirichlet(np.ones(n))  # population composition (from symmetric Dirichlet distribution)\n",
    "\n",
    "    def __repr__(self):\n",
    "        return 'Number of platform products = {:d}\\nNumber of consumer types = {:d}\\nPopulation composition = {}\\nUtilities:\\n{}'.format(self.m,self.n,self.p,self.c)\n",
    "\n",
    "    def ccp(self,u):\n",
    "        '''Conditional choice probabilities, assuming choices in row'''\n",
    "        u = np.asarray(u).reshape((self.m,self.n))   # convert to matrix\n",
    "        u = u - np.amax(u,axis=0) # de-max by column (avoid exp of large numbers)\n",
    "        e = np.exp(u)\n",
    "        esum = e.sum(axis=0)  # sums of exps\n",
    "        return e/esum         # vector of choice probabilities\n",
    "\n",
    "    def shares(self,pr):\n",
    "        '''Market shares from choice probabilities'''\n",
    "        out = pr @ self.p  # one-dim vector\n",
    "        return out[:,np.newaxis]  # column vector\n",
    "\n",
    "    def F(self,u):\n",
    "        '''Fixed point equation u=F(u)'''\n",
    "        pr = self.ccp(u)     # matrix of choice probabilities\n",
    "        sh = self.shares(pr) # market shares\n",
    "        u1 = self.c + sh     # updated utilities\n",
    "        return u1.reshape(self.m*self.n) # return one dimensional array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of platform products = 3\n",
      "Number of consumer types = 2\n",
      "Population composition = [0.07395594 0.92604406]\n",
      "Utilities:\n",
      "[[0.71361939 0.51952517]\n",
      " [0.27539451 0.30194578]\n",
      " [0.73362029 0.39378844]]\n",
      "SA from zero utilities:\n",
      "iter 0, err = 1.067e+00\n",
      "iter 1, err = 3.914e-02\n",
      "iter 2, err = 1.352e-02\n",
      "iter 3, err = 4.718e-03\n",
      "iter 4, err = 1.658e-03\n",
      "iter 5, err = 5.859e-04\n",
      "iter 6, err = 2.080e-04\n",
      "iter 7, err = 7.412e-05\n",
      "iter 8, err = 2.650e-05\n",
      "iter 9, err = 9.496e-06\n",
      "iter 10, err = 3.411e-06\n",
      "iter 11, err = 1.227e-06\n",
      "iter 12, err = 4.422e-07\n",
      "Equilibrium found!\n",
      "Equilibrium choice probabilities:\n",
      "[[0.39502849 0.39314102]\n",
      " [0.2270464  0.28174725]\n",
      " [0.37792511 0.32511173]]\n",
      "Equilibrium market shares:\n",
      "[0.39328061 0.2777018  0.32901759]\n"
     ]
    }
   ],
   "source": [
    "def printiter(**kwargs):\n",
    "    print('iter %d, err = %1.3e'%(kwargs['iter'],kwargs['err']))\n",
    "md = model(m=3,n=2)\n",
    "print(md)\n",
    "print('SA from zero utilities:')\n",
    "x = solve_sa(md.F,x0=np.zeros(md.n*md.m),callback=printiter)\n",
    "print('Equilibrium found!')\n",
    "ccp = md.ccp(x)\n",
    "shares = md.shares(ccp).squeeze()  # make one-dim array\n",
    "print('Equilibrium choice probabilities:',ccp,'Equilibrium market shares:',shares,sep='\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Further learning resources\n",
    "\n",
    "- Oscar Veliz videos on fixed point iterations\n",
    "  [https://www.youtube.com/watch?v=OLqdJMjzib8](https://www.youtube.com/watch?v=OLqdJMjzib8)\n",
    "  [https://www.youtube.com/watch?v=FyCviw2ZA2o](https://www.youtube.com/watch?v=FyCviw2ZA2o)  \n",
    "- Non-technical review of platform competition\n",
    "  [https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3502964](https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3502964)  \n",
    "- Blog posts on platform competition by Paul Belleflamme\n",
    "  [http://www.ipdigit.eu/2020/04/an-introduction-to-the-economics-of-platform-competition-part-1/](http://www.ipdigit.eu/2020/04/an-introduction-to-the-economics-of-platform-competition-part-1/)\n",
    "  [http://www.ipdigit.eu/2020/04/an-introduction-to-the-economics-of-platform-competition-part-2/](http://www.ipdigit.eu/2020/04/an-introduction-to-the-economics-of-platform-competition-part-2/)\n",
    "  [http://www.ipdigit.eu/2020/04/an-introduction-to-the-economics-of-platform-competition-part-3/](http://www.ipdigit.eu/2020/04/an-introduction-to-the-economics-of-platform-competition-part-3/)  "
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "date": 1612589585.6845,
  "download_nb": false,
  "filename": "22_succ_aprx.rst",
  "filename_with_path": "22_succ_aprx",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  },
  "title": "Foundations of Computational Economics #22"
 },
 "nbformat": 4,
 "nbformat_minor": 4
}