{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Newton's Method - Illustrative Examples\n", "\n", "This code illustrate Newton's method for various functions defined below. In particular, it contains a few examples which show the importance of every hypothesis in the theoretical quadratic convergence result." ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pylab as pl\n", "import matplotlib.pyplot as plt\n", "plt.rcParams['figure.dpi']= 100 # parameter for resolution of graphics\n", "import time\n", "\n", "v = 1 # variant corresponding to the number of the function below\n", "Niter = 10 # Number of iterations\n", "x0 = -1 # Initialization\n", "a = -2.5 # Lower bound for the plot interval\n", "b = 2.5 # Upper bound for the plot interval" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Various objective functions\n", "\n", "Here we consider multiple functions to be tested with Newton's method\n", "\n", "Case 0: $f(x) = x^2$ (quadratic function, convergence in $1$ iteration)\n", "\n", "Case 1: $f(x) = x^6/6-x^2/2+x$ (here the choice of the initialization is important)\n", "\n", "Case 2: $f(x) = x^2-\\sin x$ \n", "\n", "Case 3: $f(x) = x^2+\\exp x$\n", "\n", "Case 4: $f(x) = x^4$ (this does not verify the non-degeneracy hypothesis: quadratic convergence is not attained)\n", "\n", "Case 5: $f(x) = \\sqrt{1+x^2}$ (here the choice of initialization is important: for $|x|<1$ we have cubic convergence, while for $|x|\\geq 1$ the algorithm diverges)" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [], "source": [ "def fun(x,v): # function definition\n", " if v==0:\n", " return x**2\n", " if v==1:\n", " return x**6/6-x**2/2+x\n", " if v==2:\n", " return x**2-np.sin(x)\n", " if v==3:\n", " return x**2+np.exp(x)\n", " if v==4:\n", " return x**4\n", " if v==5:\n", " return np.sqrt(1+x**2)\n", "def der(x,v): # first derivative\n", " if v==0:\n", " return 2*x\n", " if v==1:\n", " return x**5-x+1\n", " if v==2:\n", " return 2*x-np.cos(x)\n", " if v==3: \n", " return 2*x+np.exp(x)\n", " if v==4:\n", " return 4*x**3\n", " if v==5:\n", " return x/np.sqrt(1+x**2)\n", "def der2(x,v): # second derivative\n", " if v==0: \n", " return 2\n", " if v==1:\n", " return 5*x**4-1\n", " if v==2: \n", " return 2+np.sin(x)\n", " if v==3:\n", " return 2+np.exp(x)\n", " if v==4:\n", " return 12*x**2\n", " if v==5:\n", " return 1/np.sqrt(1+x**2)**3\n", "\n", "# List of optimizers for the above functions\n", "if v==0:\n", " analytic = 0\n", "if v==1:\n", " analytic = -1.1673039782614187\n", "if v==2:\n", " analytic = 0.45018361129487355\n", "if v==3:\n", " analytic = -0.35173371124919584\n", "if v==4:\n", " analytic = 0\n", "if v==5:\n", " analytic = 0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Algorithm: Newton method in dimension one\n", "\n", "
\n", " **Initialization:** Choose the starting point $x_0$\n", "\n", " **Step $i$:** \n", " - Compute $f(x_{i-1}),f'(x_{i-1}),f''(x_{i-1})$ and approximate $f$ around $x_{i-1}$ by its second-order taylor expansion\n", "$$ p(x) = f(x_{i-1})+f'(x_{i-1})(x-x_i)+\\frac{1}{2} f''(x_{i-1})(x-x_{i-1})^2.$$\n", " - choose $x_i$ as the minimizer of the quadratic function $p$:\n", " $$ x_i = x_{i-1} - \\frac{f'(x_{i-1})}{f''(x_{i-1})}.$$\n", " - replace $i$ with $i+1$ and loop\n", "
" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAFdCAYAAACet25NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl41OW9/vH3JwESkpAQ9kASFqEgSMWFgloFEVsXrMtRrPZYpba1FhDx1AXbelrrcWkt1arVal0vq7YVixvu/hBRq4iKCLIKISEh7AlbApk8vz+SGUNIIMtMnlnu13V9r5DZcjOlzp1n+X7NOYeIiIgkriTfAURERMQvlQEREZEEpzIgIiKS4FQGREREEpzKgIiISIJTGRAREUlwKgMiIiIJTmVAREQkwakMiIiIJDiVARERkQSnMiAiIpLg2vkOUJ+ZGdAb2OE7i4iISAzqBBS7Zlx8KOrKADVFoMh3CBERkRiWC6xv6oOjsQzsACgsLCQzM9N3FhERkZhRXl5OXl4eNHN0PRrLAACZmZkqAyIiIm2gVQsIzewGM3Nmdled21LN7D4z22JmO81slpn1bH1UERERiYQWlwEzGwlcAXxe764/AWcBFwBjqFkD8FxLf46IiIhEVovKgJllAH8HfgJsq3N7FnA5cI1z7m3n3EJgEnC8mY0OQ14REZG4UF1d7TtCSEtHBu4DXnbOvVnv9mOA9kDodufcMmAdcFwLf5aIiEjcGTduHPn5+bz5Zv2P0rbX7AWEZvZ94GhgZAN39wL2Oue217u9tPa+hl4vBUipc1On5mYSERGJNStXrqS4uDgqFss3a2TAzPKAu4EfOOcqwpRhBlBW59A5BkREJK7t2bOH4uJiAAYMGOA5TfOnCY4BegCfmFmVmVVRs0jwqto/lwIdzKxzvef1BDY08pq3AVl1jtxmZhIREYkpa9asAWq20Xft2tVzmuZPE7wFDK9326PAMuAOoBDYB5wCzAIws8FAPvBBQy/onKsEKoPf15yNWEREJH6tXr0agMMOOywqPveaVQacczuAL+reZma7gC3OuS9qv38YmGlmW4Fy4B7gA+fcf8ITWUREJLZ99dVXQHRMEUBkzkA4HaimZmQgBXgN+HkEfo6IiEhMqjsyEA1aXQacc2PrfV8BTK49REREpJ5oGxlo1emIRUREpPmibWRAZUBERKQNVVdXh3YTaGRAREQkARUXF1NZWUlycjL5+fm+4wAqAyIiIm0quF6gb9++tGsXiXX8zacyICIi0oaibb0AqAyIiIi0qeDIgMqAiIhIggqODETL4kFQGRAREWlTmiYQERFJcNF2wiFQGRAREWkz5eXlbN68GVAZEBERSUjBUYFu3bqRmZnpOc3XVAZERETaSDSuFwCVARERkTYTjesFQGVARESkzWhkQEREJMFF4wmHQGVARESkzUTjCYdAZUBERKRNVFVVUVBQAGhkQEREJCGtW7eOQCBASkoKOTk5vuPsR2VARESkDdSdIkhKiq6P3+hKIyIiEqeidVshqAyIiIi0iWjdVggqAyIiIm0iWrcVgsqAiIhIm4jWbYWgMiAiIhJxzrn4GRkwsyvN7HMzK689PjCz0+vcP9fMXL3jgfDHFhERiR1btmyhvLwcgH79+vkN04B2zXx8EXADsBIw4FLgeTM7yjm3pPYxDwE31XnO7lanFBERiWErVqwAIC8vj44dO3pOc6BmjQw45150zs1xzq10zq1wzv0S2AmMrvOw3c65DXWO8rAmboEnn3ySCy+8kFmzZvmOIiIiCShYBr7xjW94TtKwFq8ZMLNkM/s+kA58UOeuH5jZZjP7wsxuM7O0Q7xOipllBg+gU0szNWbhwoX885//5P333w/3S4uIiBxSsAwMHjzYc5KGNXeaADMbTs2Hfyo1owLnOueW1t79FFAAFAPfBO4ABgPnHeQlZwD/29wczTFw4EDg65WcIiIibWn58uVA9I4MNLsMAMuBEUAWcD7wuJmNcc4tdc49WOdxi82sBHjLzA5zzjX2SXwbMLPO952oWZsQNsGVm6tWrQrny4qIiDRJ3I0MOOf2AsFP1YVmNhKYBlzRwMM/rP06EGiwDDjnKoHK4Pdm1txIhxQsA1999RXOuYj8DBERkYYEAgFWrlwJRO/IQDjOM5AEpDRy34jaryVh+Dkt1rdvX5KTk9mzZw8lJV6jiIhIgiksLKSyspIOHTrQt29f33Ea1NzzDNxmZieZWT8zG25mtwFjgb+b2WFm9mszO6b2/u8BTwDznHOfRyB7k3Xo0IH8/HxAUwUiItK2glMEAwcOJDk52XOahjV3ZKAHNR/wy4G3gJHAd51zbwB7gfHA68Ay4I/ALOCssKVtheBUgRYRiohIW4r2xYPQzDUDzrnLD3JfITCm1YkiZODAgbz55psqAyIi0qai/RwDkEDXJtCOAhER8SHadxJAApYBjQyIiEhbioVpgoQpAzrxkIiItLU9e/awbt06QGUgKgSvH71t2za2bt3qOY2IiCSCVatW4Zyjc+fOdO/e3XecRiVMGUhPT6dXr16ARgdERKRt1F08GM0nvEuYMgCaKhARkbYVC4sHIcHKgBYRiohIW4qFxYOQoGVA2wtFRKQtxMI5BiDByoCmCUREpC1pmiAKaZpARETaypYtW9iyZQvw9S+j0Sohy0BxcTG7d+/2nEZEROJZcFQgNzeX9PR0z2kOLqHKQJcuXejcuTMAX331lec0IiISz4KLB6N9igASrAyYmaYKRESkTcTK4kFIsDIA2lEgIiJtQ2UgimlHgYiItAVNE0QxTROIiEikVVdXs3LlSkAjA1FJ0wQiIhJphYWFVFZW0r59e/r16+c7ziElXBkIThMUFBSwb98+z2lERCQeBacIBg4cSHJysuc0h5ZwZSAnJ4fU1FQCgUDoGtMiIiLhtGzZMiA2pgggActAUlISAwYMADRVICIikbFkyRIAhg0b5jlJ0yRcGQDtKBARkchaunQpoDIQ1bSjQEREIsU5FxoZGDp0qOc0TZPQZUDTBCIiEm4bNmxg27ZtJCUlMWTIEN9xmqRZZcDMrjSzz82svPb4wMxOr3N/qpndZ2ZbzGynmc0ys57hj906wWkClQEREQm34BTBYYcdRmpqquc0TdPckYEi4AbgGOBY4G3geTMLTor8CTgLuAAYA/QGngtP1PAJng1q1apVVFVVeU4jIiLxJNamCKCZZcA596Jzbo5zbqVzboVz7pfATmC0mWUBlwPXOOfeds4tBCYBx5vZ6PBHb7n8/HxSU1PZu3cva9eu9R1HRETiSKztJIBWrBkws2Qz+z6QDnxAzWhBe+DN4GOcc8uAdcBxB3mdFDPLDB5Ap5ZmaqqkpKTQ3s/giSFERETCIdZ2EkALyoCZDTeznUAl8ABwrnNuKdAL2Ouc217vKaW19zVmBlBW5yhqbqaWCE4VqAyIiEi4xOJOAmjZyMByYAQwCrgfeNzMWvM3vg3IqnPktuK1mixYBoJniRIREWmtWNxJANCuuU9wzu0FgsvwF5rZSGAa8A+gg5l1rjc60BPYcJDXq6RmlAEAM2tupBYJ/o+kkQEREQmXWNxJAOE5z0ASkAIsBPYBpwTvMLPBQD41awqiiqYJREQk3GJxigCaOTJgZrcBr1CzKLATcDEwFviuc67MzB4GZprZVqAcuAf4wDn3n7CmDoPgAsLS0lK2b99O586dPScSEZFYF4uLB6H5IwM9gCeoWTfwFjCSmiLwRu3904GXgFnAPGqmB84LT9TwyszMJCcnB9DogIiIhEcsbiuEZo4MOOcuP8T9FcDk2iPqDRkyhJKSEpYvX86oUaN8xxERkRgWqzsJIEGvTRCkdQMiIhIupaWlMbmTAFQGAG0vFBGR1guOCsTaTgJI8DKg7YUiIhIusTpFAAleBoIjAytXriQQCHhOIyIisSxWdxJAgpeB/Px8UlJSdMEiERFptVjdSQAJXgaSk5N1wSIREWm1WN5JAAleBkA7CkREpPVieScBqAyoDIiISKvF8k4CUBkINThtLxQRkZaK5SkCUBnQyICIiLRaLC8eBJWBUBnYsGEDZWVlntOIiEgs+vzzzwEYPny45yQtk/BlQBcsEhGR1ggEAixevBiAESNGeE7TMglfBkBTBSIi0nKrV69m165ddOzYkUGDBvmO0yIqA6gMiIhIy3322WdAzRRBcnKy5zQtozKAdhSIiEjLLVq0CIAjjzzSc5KWUxlAIwMiItJywZGBWF0vACoDgC5YJCIiLRccGVAZiHF9+/YlNTWVyspK1qxZ4zuOiIjEiM2bN7N+/XogdrcVgsoAUHPBouBZo4InjhARETmU4KjAwIED6dSpk+c0LacyUCt41qgvvvjCcxIREYkVwfUCsbx4EFQGQo444ghAZUBERJouHtYLgMpAiMqAiIg0l0YG4kywDCxbtoy9e/d6TiMiItGusrKSL7/8EkiwkQEzm2FmC8xsh5ltNLPZZja43mPmmpmrdzwQ3tjhl5eXR6dOnaiqqmLlypW+44iISJRbunQpVVVVdOnShdzcXN9xWqW5IwNjgPuA0cCpQHvgdTNLr/e4h4CcOsd1rcwZcWamqQIREWmyulMEZuY5Tes0qww4505zzj3mnFvinFsEXAbkA8fUe+hu59yGOkd5mPJGlHYUiIhIU8XL4kFo/ZqBrNqvW+vd/gMz22xmX5jZbWaW1tgLmFmKmWUGD8DbRs3gyIDONSAiIocSL4sHAdq19IlmlgTcBbznnKv7q/RTQAFQDHwTuAMYDJzXyEvNAP63pTnCSdMEIiLSFM65uBoZMOdcy55odj9wOvBt51zRQR43DngLGOicW93A/SlASp2bOgFFZWVlZGZmtihbS5WWltKrVy/MLHRtahERkfoKCgro168f7du3Z+fOnXTo0MF3JADKy8vJysoCyGrOFH2LpgnM7F5gAnDywYpArQ9rvw5s6E7nXKVzrjx4ADtakikcevToQbdu3XDOhbaLiIiI1BecIhg6dGjUFIHWaO7WQqstAucC45xzTbmqT3D8pKS54dqadhSIiEhTxNMUATR/ZOA+4L+Bi4EdZtar9ugIYGaHmdmvzewYM+tnZt8DngDmOec+D2/0yFAZEBGRQ4mnxYPQ/AWEV9Z+nVvv9knAY8BeYDxwNZAOFAKzgFtanLCNaXuhiIgcSrAMxMvIQLPKgHPuoGdVcM4VUnNiopilkQERETmYLVu2sGZNzSx5vJQBXZugnuDIQGFhIeXlMXGuJBERaUMff/wxAIMGDSI7O9tzmvBQGagnOzubPn36ADr5kIiIHGjBggUAjBw50nOS8FEZaICmCkREpDEqAwlCZUBERBrinOOjjz4CVAbinsqAiIg0ZP369WzYsIHk5GSOOuoo33HCRmWgAdpeKCIiDQlOEQwbNoy0tEavwRdzVAYaMHToUAA2btzIxo0bPacREZFoESwD3/rWtzwnCS+VgQakp6czYMAAQDsKRETka/G4eBBUBho1fPhw4OvzT4uISGKrrq4OnWNAZSBBBBeGfPrpp56TiIhINFi1ahXbt28nNTU1tNA8XqgMNEJlQERE6gpOEYwYMYL27dt7ThNeKgONCJaBpUuXUlFR4TmNiIj4Fq/rBUBloFG5ubl07dqVQCCgLYYiIhK3OwlAZaBRZqapAhERAaCqqir0WaCRgQSjMiAiIlCzzXzPnj1kZmYyaNAg33HCTmXgIFQGREQECF2P4NhjjyUpKf4+OuPvbxRGwTKwaNEiAoGA5zQiIuJLPC8eBJWBgxo0aBBpaWns2bOH5cuX+44jIiKeqAwksOTkZI488khAUwUiIolqz549LF68GFAZSFhaNyAiktg++eQTAoEAPXv2JC8vz3eciFAZOASVARGRxDZ//nwATjjhBMzMc5rIUBk4hLplwDnnOY2IiLS19957D6gpA/FKZeAQjjjiCNq1a8e2bdtYt26d7zgiItKGqqurQ2Xg29/+tuc0kdOsMmBmM8xsgZntMLONZjbbzAbXe0yqmd1nZlvMbKeZzTKznuGN3XZSUlIYOnQooKkCEZFEs2zZMrZu3UrHjh1DI8XxqLkjA2OA+4DRwKlAe+B1M0uv85g/AWcBF9Q+vjfwXOuj+nP00UcDKgMiIokmOCowatSouLtSYV3NKgPOudOcc48555Y45xYBlwH5wDEAZpYFXA5c45x72zm3EJgEHG9mo8Mbve1oEaGISGIKLh6M5ykCaP2agazar1trvx5DzWjBm8EHOOeWAeuA41r5s7xRGRARSUyJsHgQWlEGzCwJuAt4zzkXvMZvL2Cvc257vYeX1t7X0OukmFlm8AA6tTRTpARPPFRUVMTmzZs9pxERkbZQUlLC6tWrMTOOOy5mf59tktaMDNwHHAF8v5UZZgBldY6iVr5e2GVmZjJw4EBAowMiIokiOCowfPhwsrKyDvHo2NaiMmBm9wITgJOdc3U/vDcAHcysc72n9Ky9ryG3UTPdEDxyW5Ip0jRVICKSWBJhS2FQc7cWWm0ROBcY55xbU+8hC4F9wCl1njOYmkWGHzT0ms65SudcefAAdjQnU1sJ7ihYuHCh5yQiItIWEmXxIEC7Zj7+PuBi4Gxgh5kF1wGUOef2OOfKzOxhYKaZbQXKgXuAD5xz/wlbag+CF6f48MMPPScREZFI27lzZ2gkON4XD0LzpwmupGYofy5QUue4sM5jpgMvAbOAedRMD5zX2qC+jRw5EjOjoKCA0tJS33FERCSCPvroIwKBAHl5eeTn5/uOE3HNPc+ANXI8VucxFc65yc65Ls65dOfcec65xtYLxIzMzEwOP/xwQKMDIiLxLpGmCEDXJmiWUaNGASoDIiLxLlHOLxCkMtAMKgMiIvGvqqqK999/H9DIgDQgWAYWLFhAdXW15zQiIhIJixcvZufOnWRmZnLEEUf4jtMmVAaa4YgjjiAtLY3y8nKWLVvmO46IiERAcL3AcccdR3Jysuc0bUNloBnatWvHMcccA2iqQEQkXr399tsAjBkzxnOStqMy0ExaNyAiEr8CgQBz584F4JRTTjn4g+OIykAzqQyIiMSvTz75hO3bt5OVlRUaCU4EKgPNFCwDixcvZvfu3Z7TiIhIOL311lsAjB07NmHWC4DKQLPl5uaSk5NDIBDQdQpEROJMsAwk0hQBqAw0m5kxevRoQFMFIiLxpKKiIrSTQGVADknrBkRE4s8HH3xARUUFOTk5odPPJwqVgRZQGRARiT/BKYJx48ZhZp7TtC2VgRY49thjSUpKorCwkJKSEt9xREQkDBJ1vQCoDLRIRkYGw4YNAzQ6ICISD8rLy1mwYAGgMiDNoKkCEZH48c477xAIBBg4cCD5+fm+47Q5lYEWUhkQEYkfiTxFACoDLRbcXvjRRx+xb98+z2lERKQ1VAakRYYOHUp2dja7du3i008/9R1HRERaqLS0lC+++AKAk08+2XMaP1QGWigpKYkTTzwRgHnz5nlOIyIiLRW8SuGIESPo1q2b5zR+qAy0wkknnQSoDIiIxLK65xdIVCoDrRAsA++++y7V1dWe04iISHM553j11VcBOPXUUz2n8UdloBWOOuoo0tPT2b59e2i+SUREYseiRYtYv349aWlpjB071nccb1QGWqFdu3accMIJgKYKRERi0csvvwzA+PHjSU1N9ZzGH5WBVtK6ARGR2PXSSy8BcOaZZ3pO4lezy4CZnWRmL5pZsZk5Mzun3v2P1d5e93g1fJGjS90y4JzznEZERJpq06ZNoRPHnXHGGZ7T+NWSkYF0YBEw+SCPeRXIqXNc1IKfExNGjhxJSkoKpaWlrFy50nccERFpoldeeQXnHCNGjCA3N9d3HK+aXQacc684537lnPv3QR5W6ZzbUOfY1oqMUS01NTV0amJNFYiIxI7geoFEnyKAyK0ZGGtmG81suZndb2ZdG3ugmaWYWWbwADpFKFPEaN2AiEhs2bdvH6+99hoAEyZM8JzGv0iUgVeBHwKnANcDY4BXzCy5kcfPAMrqHEURyBRRKgMiIrHlvffeo6ysjG7dujFy5EjfcbwLexlwzj3jnHvBObfYOTcbmACMBMY28pTbgKw6R8xN3Bx33HEkJydTUFBAQUGB7zgiInIIwSmCM844g+Tkxn5XTRwR31ronPsK2AwMbOT+SudcefAAdkQ6U7hlZGRw9NFHAzVnIxQRkeimLYX7i3gZMLNcoCtQEumf5ZOmCkREYsPq1atZtmwZycnJfOc73/EdJyq05DwDGWY2wsxG1N7Uv/b7/Nr7/mBmo82sn5mdAjwPrAJeC2fwaKMyICISG4JTBCeeeCKdO3f2nCY6tGRk4Fjg09oDYGbtn28GAsA3gReAFcDDwELgROdcZavTRrFvf/vbACxfvpwNGzZ4TiMiIo3RlsIDteQ8A3Odc9bAcZlzbo9z7rvOuR7OuQ7OuX7OuZ8650ojET6adOnShREjagZLgpfDFBGR6FJWVsbcuXMBlYG6dG2CMArOPb3xxhuek4iISENefPFF9u7dy5AhQxgyZIjvOFFDZSCMgmXg9ddf13UKRESi0L/+9S8AJk6ciJl5ThM9VAbC6IQTTqBjx46UlJSwZMkS33FERKSOsrIyXn215rp5F1xwgec00UVlIIxSU1MZM2YMUDM6ICIi0aPuFMGwYcN8x4kqKgNhVneqQEREooemCBqnMhBmwTLwzjvvUFFR4TmNiIiApggORWUgzIYOHUrv3r2pqKhg/vz5vuOIiAhfTxEcfvjhmiJogMpAmJmZpgpERKJMcIrgggsu0BRBA1QGIkBlQEQkemiK4NBUBiJg/PjxACxatEinJhYR8UxTBIemMhAB3bt3D13S+M033/ScRkQksWmK4NBUBiJEUwUiIv5piqBpVAYiRKcmFhHx77nnntMUQROoDETI8ccfT1paGqWlpSxevNh3HBGRhPToo48C8MMf/lBTBAehMhAhKSkpjB07FoBXXnnFbxgRkQS0atUq3n33XZKSkrjkkkt8x4lqKgMRNGHCBACef/55z0lERBLP448/DtRM2/bp08dzmuimMhBB3/ve9wD4z3/+Q2lpqec0IiKJIxAIhMrApEmTPKeJfioDEdSnTx+OPfZYnHO8+OKLvuOIiCSMt99+m8LCQrKzs0O/mEnjVAYi7JxzzgFg9uzZnpOIiCSO4MLBiy66iNTUVM9pop/KQISdffbZQM3Jh3bu3Ok5jYhI/Nu+fTv//ve/AU0RNJXKQIQNGzaMAQMGUFlZqRMQiYi0gX/84x9UVFRwxBFHcMwxx/iOExNUBiLMzEKjA9pVICISecEpgkmTJuncAk2kMtAGgmXgpZdeoqqqynMaEZH49eWXX/Lhhx+SnJzMD37wA99xYkazy4CZnWRmL5pZsZk5Mzun3v1mZjebWYmZ7TGzN81sUPgix54TTjiBrl27snXrVubPn+87johI3Prb3/4GwJlnnknPnj09p4kdLRkZSAcWAZMbuf864CrgZ8AoYBfwmpkl7HLOdu3a6QREIiIRtnPnTh5++GEAfvrTn3pOE1uaXQacc684537lnPt3/fusZnLmauAW59zzzrnPgR8CvYFz6j8+kdRdN6ALF4mIhN+TTz5JWVkZAwcO5PTTT/cdJ6aEe81Af6AX8GbwBudcGfAhcFxDTzCzFDPLDB5ApzBnigrf+c53SE1NZc2aNbpwkYhImDnn+POf/wzA1KlTSUrSkrjmCPe71av2a/1z75bWua++GUBZnaMozJmiQnp6OuPHjwc0VSAiEm5vvfUWX375JRkZGVx22WW+48ScaKhOtwFZdY5cv3EiJ3g2wmeffdZzEhGR+HL33XcDNdsJMzMzPaeJPeEuAxtqv9Zfwtmzzn37cc5VOufKgwewI8yZosZ5551H+/bt+fzzz/niiy98xxERiQurVq3i5ZdfBmDKlCme08SmcJeBNdR86J8SvKF2HcAo4IMw/6yYk52dzRlnnAHAU0895TmNiEh8uO+++3DOcfrpp/ONb3zDd5yY1JLzDGSY2QgzG1F7U//a7/NdzTL5u4Bfmdn3zGw48ARQDOhKPcDFF18M1JQB7SoQEWmdHTt28MgjjwBw1VVXeU4Tu1oyMnAs8GntATCz9s83137/e+Ae4EFgAZABnOacq2hd1PgwYcIEMjIyKCgo4IMPEn6wRESkVZ544gnKy8v5xje+wXe+8x3fcWJWS84zMNc5Zw0cl9Xe75xzNznnejnnUp1z451zK8KePEalpaVx7rnnApoqEBFpjaqqKu666y5A2wlbS++cB8HzZf/zn/9k3759ntOIiMSmZ555hlWrVtG1a1dtJ2wllQEPTjnlFLp3786mTZt48803D/0EERHZTyAQ4JZbbgHgf/7nf8jIyPCcKLapDHjQrl07LrzwQkBTBSIiLfHss8+yfPlysrOzmTy5sUvlSFOpDHgS3FXw73//m927d3tOIyISO6qrq/nd734HwNVXX62TDIWByoAno0ePpn///uzatYsXX3zRdxwRkZgxe/ZslixZQmZmprYThonKgCdmtt85B0RE5NCcc6FRgauuuorOnTt7ThQfVAY8CpaBOXPmsGFDg2drFhGROl566SU+++wzMjIyuPrqq33HiRsqAx4NHTqU0aNHU1VVxaOPPuo7johIVHPOcfPNNee3mzx5Ml27dvWcKH6oDHh2xRVXAPDQQw9RXV3tOY2ISPR69tln+fjjj0lLS+Oaa67xHSeuqAx4NnHiRLKyslizZg1vvPGG7zgiIlGpoqKC66+/HoDrrruOHj16eE4UX1QGPEtLS+OHP/whAH/96189pxERiU733HMPa9asoXfv3vziF7/wHSfuqAxEgeBUwQsvvEBxcbHnNCIi0WXTpk2hsw3eeuutpKene04Uf1QGosCwYcM44YQTCAQCoUtxiohIjd/85jeUl5dz9NFHc8kll/iOE5dUBqJE3YWEgUDAcxoRkeiwdOnS0BTqH//4R12ZMEL0rkaJ888/n+zsbNatW8drr73mO46ISFS49tprCQQCnHPOOYwdO9Z3nLilMhAlOnbsyKWXXgrAAw884DmNiIh/c+bMYc6cObRr147f//73vuPENZWBKPLTn/4UgJdffpl169Z5TiMi4s+OHTvDsjkjAAAV4ElEQVT42c9+BsC0adMYNGiQ50TxTWUgihx++OGcfPLJVFdXc/fdd/uOIyLizYwZMygsLKR///789re/9R0n7qkMRJlrr70WgAcffJBt27Z5TiMi0vbmz5/PX/7yF6Dmv4XaShh5KgNR5rTTTmP48OHs3LlTawdEJOFUVFTw4x//GOccP/rRjxg/frzvSAlBZSDKmBnXXXcdAHfffTcVFRWeE4mItJ3f/e53LF++nF69enHnnXf6jpMwVAai0IUXXkheXh6lpaU88cQTvuOIiLSJzz77LLRr4C9/+QvZ2dmeEyUOlYEo1L59+9AVue68806dhEhE4t6uXbv4wQ9+QFVVFf/1X//Fueee6ztSQgl7GTCz35iZq3csC/fPiXc//vGPyc7OZuXKlcyePdt3HBGRiJoyZQpLly4lJycntHhQ2k6kRgaWADl1jm9H6OfErYyMDCZPngzAHXfcgXPOcyIRkch4/PHHeeyxx0hKSuKpp57S5Yk9iFQZqHLObahzbI7Qz4lrU6dOJTU1lQULFvDOO+/4jiMiEnZLly7l5z//OVBzQSKdctiPSJWBQWZWbGZfmdnfzSy/sQeaWYqZZQYPoFOEMsWcHj16MGnSJAB+/etfa3RAROLK7t27mThxIrt37+aUU07hxhtv9B0pYUWiDHwIXAacBlwJ9AfeNbPGPuRnAGV1jqIIZIpZN954I6mpqcyfP585c+b4jiMiEhbOOSZPnsySJUvo1asXf//730lOTvYdK2GFvQw4515xzv3LOfe5c+414AygMzCxkafcBmTVOXLDnSmW5ebmMnXqVKDm9JzV1dWeE4mItN4dd9wRWifw97//nZ49e/qOlNAivrXQObcdWAEMbOT+SudcefAAdkQ6U6y54YYbyMrKYvHixTz11FO+44iItMozzzzDjBkzgJqTq40bN85zIol4GTCzDOAwoCTSPytedenSheuvvx6oWTuwd+9ez4lERFpm/vz5ocu1T58+nSlTpnhOJBCZ8wzcaWZjzKyfmR0P/BsIAE+H+2clkmnTppGTk8PatWv561//6juOiEizrVixgrPPPpu9e/dy7rnn8oc//MF3JKkViZGBXGo++JcD/wS2AKOdc5si8LMSRlpaGjfddBNQc+7uHTs0myIisWP9+vWcccYZbN26lW9961s8+eSTWjAYRSKxgPD7zrnezrkU51xu7ferw/1zEtHll1/OwIED2bRpEzNnzvQdR0SkSYqKihg7diyrV6+mf//+vPDCC6SlpfmOJXXo2gQxpH379vzf//0fULMSd82aNZ4TiYgcXGFhIWPHjmXVqlX069ePt99+WzsHopDKQIy54IILOPnkk9mzZw+TJ0/WiYhEJGoVFBQwZswYVq9ezYABA3jnnXfo16+f71jSAJWBGGNm3H///XTo0IFXXnmFWbNm+Y4kInKAVatWMXbsWNasWcNhhx3G3Llzyc9v9GS04pnKQAwaPHhwaKvhtGnTKC8v95xIRORr8+bNY9SoUaxdu5ZBgwYxd+5c8vLyfMeSg1AZiFE33ngjAwcOpLi4mF//+te+44iIAPDYY48xfvx4tm7dysiRI3nnnXfIzdWJZaOdykCMSk1NDV3z+95772XhwoWeE4lIIquurmbGjBlMmjSJffv2cf755zN37lxycnJ8R5MmUBmIYaeeeioXXXQR1dXVXHHFFezbt893JBFJQBs3buSss87i9ttvB+BXv/oV//jHP7R9MIaoDMS4mTNn0rlzZxYuXMhvf/tb33FEJMG8+uqrfPOb32TOnDmkpKTwxBNP8Lvf/Y6kJH28xBL9rxXjevXqxYMPPgjArbfeyty5c/0GEpGEUFlZyTXXXMPpp59OaWkpw4YNY8GCBVxyySW+o0kLqAzEgQsuuIDLL78c5xz//d//zZYtW3xHEpE49u6773Lsscfypz/9CYApU6awYMEChg8f7jmZtJTKQJy4++67GTx4MOvXr+fHP/6xTkYkImFXWlrKpZdeykknncQXX3xBt27dePHFF7nnnnvo2LGj73jSCioDcSI9PZ2nn36aDh06MHv2bF3ZUETCZt++fdx3330MHjyYJ554AoCf/OQnLFu2jAkTJnhOJ+GgMhBHjjrqqNBq3unTp/PZZ595TiQisayqqopHH32UwYMHM2XKFMrKyjjqqKP44IMPePDBB+natavviBImKgNxZtq0aZxxxhlUVFQwYcIEiouLfUcSkRhTVVXFk08+yeGHH86PfvQj1qxZQ48ePbj33ntZsGABo0eP9h1RwsyibW7ZzDKBsrKyMjIzM33HiUnbt2/n+OOP58svv+Too49m3rx5pKen+44lIlFu06ZNPPTQQ9x///0UFRUB0K1bN66//nquvPJK/XckBpSXl5OVlQWQ5Zxr8rnqVQbi1FdffcWoUaPYvHkz55xzDrNmzdK+XxE5gHMuNOz/zDPPUFlZCUD37t2ZPn06U6ZMoVOnTp5TSlOpDMgB3n//fcaNG0dlZSXXXnstv//9731HEpEosWTJEp566imeeuop1q5dG7r92GOPZerUqUycOJHU1FR/AaVFVAakQU8//TQXX3wxAH/+85+ZOnWq50Qi4kNVVRXvv/8+c+bM4aWXXmLJkiWh+zIyMjjvvPO48sorGTVqFGbmMam0RkvLQLvIRZJocNFFF7Fq1SpuuukmrrrqKgAVApEEUF1dzZIlS5g3bx7vvPMOb7zxBtu3bw/d3759e04//XQuvvhizjrrLF1HIMGpDCSAX/3qV+zatYs77riDq666iurqaqZNm+Y7loiEiXOOkpISPvnkEz755BMWLFjAe++9x7Zt2/Z7XNeuXTnttNM488wzOe2008jOzvaUWKKNykACMDNuu+02kpOTufXWW7n66quprq5m+vTpvqOJSDNUV1dTUlLCihUr+PLLL1m2bBlffvklixcvprS09IDHp6enc9xxx3HiiScyfvx4Ro0aRXJysofkEu1UBhKEmXHLLbeQlJTELbfcwjXXXENFRQU33HCD5gdFokAgEGDLli2UlpZSXFzM+vXrKS4upqioiLVr17JmzRoKCgpCq/3rS0pKYujQoRx11FEcffTRnHDCCYwYMYL27du38d9EYpHKQAIxM26++WaSkpK4+eabufHGG1mxYgUPPPAAKSkpvuOJxLRAIMCuXbvYuXMnO3bsCH0tKyvb79i+fTtbt25l69atbNu2jc2bN7Np0yY2b97cpGuKJCcn069fPw4//PDQMXToUIYPH655f2mxiO0mMLPJwLVAL2ARMNU591ETnqfdBG3gnnvuCU0XHH/88Tz33HP07NnTdyyRiKiqqmLXrl3s3r37kMeePXsO+n3w2LVr135HRUVFWLJ27dqV3r1706dPn9DRt29f+vfvT//+/cnNzaVdO/0eJw2Lqq2FZnYh8ATwM+BD4GrgAmCwc27jIZ6rMtBGXn/9dSZOnEhZWRn5+fm88MILHHnkkb5jSYJzzlFRUUF5eTllZWWUl5dTXl7Ojh07Ql/rHjt37jzgqPshvXPnTvbt29dm+ZOSkujUqROdOnUiIyODrKwssrKyyMzMJCsri+zsbLp06RL62qVLF3r27En37t3p1q2bPuilVaKtDHwILHDOTan9PgkoBO5xzt1+iOeqDLSh5cuX873vfY8VK1bQsWNHbr/9dqZMmaKzFUqrVFdXU1ZWFhoODw6JB79u376dbdu2hf5cdwi9rKwsYh/eZkZ6ejppaWkHHB07diQ9PZ2OHTse9Pu6z0lPTw8dGRkZpKenk5qaqnU44k3UlAEz6wDsBs53zs2uc/vjQGfn3NmHeL7KQBvbtm0bF198Ma+++ioA48aN45FHHqFv376ek0k0cM6xc+dONm7cyKZNm/Y7tmzZwubNm9m8eXPoz1u2bGHbtm1UV1e36ueaGZ06dQr9Rh38c/C37uBv3p06dSI9PT30NfihXPdDOvjBnZKSog9qiWvRdNKhbkAyUH+fSykwpP6DzSwFqLt6TSfBbmPZ2dnMmTOHBx54gF/84he8/fbbDB8+nLvuuotJkybpP55xyDnH9u3b2bBhw35HaWkppaWlbNy4MfR106ZNLZ4PT09PDw2FZ2dnh4bGO3fuHPo+KyuLzp0707lz59CQelZWFhkZGRqhEmkjkRgZ6A2sB453zn1Q5/bfA2Occ6PqPf43wP/Wfx2NDPixatUqLrvsMt577z0AjjvuOO68806OP/54z8mkKZxzbN26NbQtre5RUlISOjZs2MDevXub9dppaWn06NGDbt260b17d7p3707Xrl1DX7t160bXrl1DR5cuXbRLRaSNxew0QSMjA0UqA/4EAgFmzpzJb37zG3bv3g3A+eefz+23385hhx3mOV3iqq6uZsOGDRQWFlJUVERRURHr16/f78/r169vdB96Qzp37kyvXr3o1asXPXv2bPAIfvDr8rUi0S9qygCEFhB+5JybWvt9ErAOuFcLCGNHcXExN910E4888gjOOdq3b8+ll17K9OnTGTp0qO94caesrIyCggLWrVtHYWEh69atCx2FhYWsX7+eqqqqJr1W165d6dOnD71796Z3797k5OSEvgaPnj176qp0InEm2srAhcDjwBXAR9RsLZwIDHHOHXjOzP2fqzIQZT7//HOuu+46XnvttdBtp512Gtdccw3jx4/XmoImcM6xadMm1q5dy9q1aykoKKCgoIC1a9eybt06CgoKKC8/9P9vk5KS6N27N3l5eeTm5oaOunvSe/fureF5kQQVVWUAwMym8PVJhz4DrnLOfdiE56kMRKn58+czc+ZMZs+eHTpT2qBBg7jooou46KKLGDLkgPWhTRbYG+DdvyymZPVucg5L48SfDye5QzKBvQGeu+ENiv/0MmlUUGEd+f6iG+k+vNehXzMA774LJSWQkwMnngiROi17dXU1paWloQ/4ukfwg3/Pnj2HfJ2uXbuSl5dH3759yc/PJy8vj7y8vNCfc3JytA9dRBoVdWWgpVQGot/q1av585//zMMPP8yuXbtCt48YMYKJEydy6qmnctRRRzX5gijPXfcfps3MpyjQO3RbbnIxk4/+D1kLXuNsXqQ3JaH7isnhBc7iZ+6vjb/mczBtGhQVfX1bbi7cfTecd14z/rK1AoEAxcXF+/1GX/frwc4ZH2Rm9O7dm379+tG3b98Djry8PDIyMpofTkSklsqAtLkdO3bw/PPP8/TTT/P666/vN5+dnZ3NySefzLhx4xg5ciTDhw+nY8eOB7zGc9f9h/P/8C1q/hV+vY1sFO9zKY9zBQ8BjrobzKqpmZZ4kJ80WAieew7OPx/q/9MOzmY8++z+hcA5x5YtWygqKqKwsHC/IzhnX1RURCAQOOj7kZSURJ8+fejfvz99+/YNfegHv8/Ly6NDhw4HfQ0RkdZQGRCvtmzZwqxZs3j55ZeZO3fuAfPfycnJDBkyhBEjRjB48GAGDBhAfp98Lj7lMIqqe1G3CBgB+rOId5lAL0poaKd5NUYJOXT4fOF+UwaBAPTrB0VFDmhoLYMjI2M7p532M0pK1oe24DVlm127du1CQ/h1f7vv168f/fr1Izc3V1eIExGvVAYkalRVVfHxxx/z1ltvMW/ePD799FM2bdp0wOOOZAyLmNvA7Z9xKY8wnXsO+bPuZjJvnLmWXbt2UVZWxoYNQygpeaoJKccC7+x3S48ePUIL8+rO1Qfn73v16qVrwYtIVIumMxBKgmvXrh2jR49m9OjR/PKXv8Q5R3FxMZ9++imLFi1i9erVrFmzhpSFg1i048Dnp7GbdA692A6gIxW8/PLLdW4Z3KTnXXLJ9Zx11uTQ6vucnBytwBeRhKUyIBFnZqFtbxMmTAjdPveuz3ht+oGP300auzhwfUFD9pDK3/72N9LT08nKymL16jymTj308370o9MZO7aJfwERkTinaQLxJrA3QL+0UtYHeu23RHD/NQMbSOLAf6OHWjOwfv2BCwihZhFhbi6sWRO5bYYiIr60dJpAVwERb5I7JHP3NesAML6+wp0jme5U8DxnAV/vHggKfv8iEw4430Bycs32Qfh690BQ8Pu77lIREBGpS2VAvDrv96N59tqP6JO8Yb/bi5P7Ujbyu/yVn7CB/T/wS8hpdFsh1GwbfPZZ6NNn/9tzcw/cVigiIpomkCgR62cgFBGJBtpaKCIikuC0ZkBERERaRGVAREQkwakMiIiIJLioPelQU67tLiIiIl9r6WdnNC4g7AMUHfKBIiIi0phc59z6pj44GsuAAb2BBs5anxA6UVOGcknc9yCc9H6Gn97T8NL7GX6J/p52AopdMz7go26aoDZ8k9tMvLGvT5u3oznbQqRhej/DT+9peOn9DD+9pzT776wFhCIiIglOZUBERCTBqQxEn0rgt7VfpfX0foaf3tPw0vsZfnpPmynqFhCKiIhI29LIgIiISIJTGRAREUlwKgMiIiIJTmVAREQkwakMRCkz62dmD5vZGjPbY2arzey3ZtbBd7ZYZma/NLP3zWy3mW33nSfWmNlkM1trZhVm9qGZfct3plhmZieZ2YtmVmxmzszO8Z0plpnZDDNbYGY7zGyjmc02s8G+c8UClYHoNYSa/32uAIYB04GfAbf6DBUHOgD/Au73HSTWmNmFwExqtmwdDSwCXjOzHl6DxbZ0at7Hyb6DxIkxwH3AaOBUoD3wupmle00VA7S1MIaY2bXAlc65Ab6zxDozuwy4yznX2XeWWGFmHwILnHNTar9PAgqBe5xzt3sNFwfMzAHnOudm+84SL8ysO7ARGOOcm+c7TzTTyEBsyQK2+g4hiad2euoY4M3gbc656trvj/OVS+QQsmq/6r+bh6AyECPMbCAwFfir7yySkLoByUBpvdtLgV5tH0fk4GpHru4C3nPOfeE7T7RTGWhjZnZ77UKhgx1D6j2nD/Aq8C/n3EN+kkevlrynIhL37gOOAL7vO0gsiLpLGCeAPwKPHeIxXwX/YGa9gf8HvA/8NHKxYlqz3lNpkc1AAOhZ7/aewIa2jyPSODO7F5gAnOScK/KdJxaoDLQx59wmYFNTHls7IvD/gIXApNo5WqmnOe+ptIxzbq+ZLQROAWZDaBj2FOBen9lEgszMgHuAc4Gxzrk1niPFDJWBKFVbBOYCBcAvgO41/87BOaffxFrIzPKBLkA+kGxmI2rvWuWc2+kvWUyYCTxuZh8DHwFXU7M17lGvqWKYmWUAA+vc1L/23+RW59w6T7Fi2X3AxcDZwA4zC65nKXPO7fEXK/ppa2GUqt361uB/ZJ1z1rZp4oeZPQZc2sBdJzvn5rZtmthjZlOAa6lZNPgZcJVz7kO/qWKXmY2lZvSvvsedc5e1bZrYV7s9syGTnHOPtWWWWKMyICIikuC0m0BERCTBqQyIiIgkOJUBERGRBKcyICIikuBUBkRERBKcyoCIiEiCUxkQERFJcCoDIiIiCU5lQEREJMGpDIiIiCQ4lQEREZEEpzIgIiKS4P4/XQshwKBp0RIAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "uplim = max(fun(a,v),fun(b,v))+1 # set limits for the plot window\n", "dnlim = -3\n", "t1 = np.linspace(a,b,100) # Create a discretization to be used with the plots\n", "plt.figure(1)\n", "plt.ylim([dnlim,uplim]) # Set upper bounds for the figure\n", "\n", "plt.plot(t1,fun(t1,v),'k') # Plot the function to be optimized on the interval [a,b]\n", "vals = np.array([x0]); # Create an array which holds the optimization history\n", "\n", "\n", "for i in range(0,Niter):\n", " # Define the Newton interpolating polynomial: see the course \n", " def interp(x):\n", " return fun(x0,v)+der(x0,v)*(x-x0)+0.5*der2(x0,v)*(x-x0)**2\n", "\n", " plt.plot(x0,fun(x0,v),'bo')\n", " #plt.plot(t1,interp(t1),'r') # Uncomment this if you want to see the interpolating function\n", "\n", " newx = x0-der(x0,v)/der2(x0,v) # Compute the next point\n", " plt.plot(newx,fun(newx,v),'ro') # Plot the next point\n", " #plt.rc(\"savefig\", dpi=300) # Uncomment if you want to save a picture\n", " x0 = newx \n", " actval = np.array([x0]) \n", " vals = np.append(vals,actval) # Update the list of values" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After the optimization loop, the convergence history and the order of convergence is computed. For each one of the functions a value $x^*$ close to the analytical solution is provided. The order of convergence is obtained by plotting the ratio of two consecutive errors $\\displaystyle\\frac{|x_{n+1}-x^*|}{|x_n-x^*|}$ is plotted in log-log scale." ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Solution Found: -1.1673039782614187\n", "History:\n", "0 -1.0\n", "1 -1.25\n", "2 -1.1784593935169048\n", "3 -1.16753738939611\n", "4 -1.1673040828230083\n", "5 -1.1673039782614396\n", "6 -1.1673039782614187\n", "7 -1.1673039782614187\n", "8 -1.1673039782614187\n", "9 -1.1673039782614187\n", "10 -1.1673039782614187\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhkAAAFjCAYAAACQSxmxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XlclVXix/HPw77vi4DARSQVES1Nc0HcuGqajpPWT8vJNZwWbZosNSuzktJqLM1xTZnJltHUNEsvpqSYuSsqiriwuACy7+t9fn88chVXQC4X5LxfL1+9zvM89znncgO+nHOecyRZlhEEQRAEQahvRoZugCAIgiAIDycRMgRBEARB0AsRMgRBEARB0AsRMgRBEARB0AsRMgRBEARB0AsRMgRBEARB0AsRMgRBEARB0AsRMgRBEARB0AsRMgRBEARB0AsRMgRBEARB0AsRMgRBEARB0IsmFzIkSdooSVK2JEnrDd0WQRAEQRDursmFDOAL4G+GboQgCIIgCPdmYugG1JYsy9GSJPWp7eskSZIATyC/3hslCIIgCA8/W+CKXIvt2xs0ZEiS1BuYDnQGPIARsixvuuWal69f0wI4Drwqy/KBeqjeE7hUD/cRBEEQhOaqJXC5phc3dE+GNUpw+BrYcOtJSZKeBT4HpgD7gdeA7ZIktZFlOf0B684HSElJwc7O7gFvJQiCIAjNR15eHt7e3lDL0YAGDRmyLP8K/AqgjF7c5nVghSzLq69fMwUYAkwAPq5NXZIkmQPmNx2yBbCzsxMhQxAEQRAaQKOZ+ClJkhnKMMqOqmOyLGuvl7vX4ZYzgdyb/omhEkEQBEFoQI0mZAAugDGQdsvxNJT5GQBIkrQDWAc8KUnSJUmS7hZAIgD7m/61rPcWC4IgCIJwV03x6ZIBNbyuFCitKt9leEYQBEEQBD1pTD0ZGUAl4H7LcXcgteGbIwiCIAjCg2g0PRmyLJdJknQY6A9sApAkyeh6eXFDt6eiooKysrKGrlZoAGZmZpiYNJr/9QVBEB5aDb1Ohg3Q+qZDfpIkdQKyZFlORnl8NVKSpEPAAZRHWK2B1Q3VRlmWSU5OJiMjo6GqFAzAxcUFHx8fMYwmCIKgRw3951wXYNdN5c+v/zcSGCfL8g+SJLkCc1Emex4DBsmyfOtkUL2pChheXl7Y2NhgZNSYRpSEB6XVaikoKODyZWUtGV9fXwO3SBAE4eHV0OtkRAP3/NNRluXFGGB4BJQhkqqA0aJFi/u/QGiSbGxsALh8+TJ5eXkEBgZibGxs4FYJgiA8fMSf6TepmoNR9UtIeHhVfcZ79uwhOjqaWizFLwiC0Dg0gZ9bImTcgRgiefhVfcbW1tacOHFCzMERBKHpKMmDzVPht7mGbsl9id+mQrNmbW1NcXEx6ekPujWOIAhCAzi3A5Z0hyOR8MeXkNu4F7MWz/EJzVrV0yUVFRUGbokgCMI9FOeA5m04+o1SdlTBsMVg37gXsxYhQ2j2xGOsgiA0amc1sGUa5F8BJOgWDv3fBTNrQ7fsvsRwSTOyZs0aJEm6678///zT0E0UBEEQqhRnw8a/w7ejlIDh1ArG/wKDP2kSAQNET0azNHfuXPz8/G473rp16ztcLQiCIDS4+F9hy2tQkApI0P1l6Ps2mFkZumW1IkJGMzR48GC6dOlS4+srKirQarWYmZnddq6kpAQzM7MHeiKnPu4hCILwUCjKgm0zIPYHpezcGoYvAZ9uhm1XHYmf6kI1iYmJSJLEp59+ysKFC/H398fc3Jy4uDiio6ORJInvv/+e2bNn4+XlhZWVFXl5eQBcuHCBUaNG4eTkhJWVFU888QRbt26tdv973aO8vJz333+fgIAALCwscHZ2plevXkRFRRniSyEIgtCwTv8MX3VTAoZkBD2mwpSYJhswQPRkNEu5ubm3rQshSRLOzs668urVqykpKeHFF1/E3NwcJycncnJyAPjggw8wMzPjjTfeoLS0FDMzM9LS0ujRowdFRUVMnToVZ2dnIiMjGTZsGOvXr2fEiBHV6rvTPebMmUNERASTJk2ia9eu5OXlcejQIY4cOUJYWJj+vzCCIAiGUJgJv74JJ9crZZc28Jcl0LLmPc6NlQgZ9yHLMsXllYZuRjWWpsYP9ETEgAEDbjtmbm5OSUmJrnzp0iXOnTuHq6ur7tiFCxcAZXjj0KFDWFpa6s69/fbbpKWlsWfPHnr16gXA5MmTCQ4O5vXXX2f48OHVhkPudI+tW7fy5JNPsnz58jq/N0EQhCYl7ifY+k8ovKb0XvScBqEzwNTC0C2rFyJk3EdxeSWB7243dDOqiZs7ECuzun90X331FY888ki1Y7fu3fH0009XCxg3e+GFF6qFA4BffvmFrl276gIGKEt3v/jii8ycOZO4uDiCgoLueQ8HBwdOnTpFQkICAQEBdXpvgiAITULBNfjlDYjbpJRd28FfvgKvzjV6ubaoCNncDGPjxv1rvHG3TtCLrl273nfi552ePrnXuaSkJLp1u33csF27drrzN4eMO91j7ty5DB8+nEceeYSgoCAGDRrE2LFjCQ4OvmdbBUEQmgxZhlMblYBRlAmSMfT6B4S+CSbm9315TkkO0Wd+wW7GQjJa2jFquQbplknzWq3MtYJS3O0M3xsiQsZ9WJoaEzd3oKGbUY2lqf53DL21l6Gm5x7k/r179+b8+fP89NNPaDQaVq5cyb/+9S+WLl3KpEmTHrhOQRAEgypIh62vw+ktStmtvdJ74fnoPV+WXZLNb8m/oUnUcDJxPzO+K8PjKtik5ZOdnICTqo3u2opKLTM2nCAmIYN1U7rj7WTYR15FyLgPSZIeaGiiufD19SU+Pv6242fOnNGdrwknJyfGjx/P+PHjKSgooHfv3syZM0eEDEEQmi5ZhpM/Kr0XxdlgZAIhb0DIP8Hk9qUBALJKsnTB4mDqQSrlSqyLZd7+vpLWqVBua4HLks+qBYzSikpe+/4Yv55MxdhI4uTlXBEyhIfDk08+ycKFC9m3bx/du3cHoLCwkOXLl6NSqQgMDLzvPTIzM6s94WJjY0Pr1q1JSUnRW7sFQRD0Kj8Vfn4d4q8/zt+ig7Luhcftw8CZxZk3gkXaQbSyVnfuUYtHePWHTGxS0zB2dMRvzWos2twIGMVllYR/c5jdZ69hZmzEl6MfZVBQC72/vfsRIaMZ+vXXX3U9DDfr0aNHnRfEmjFjBt999x2DBw9m6tSpODk5ERkZycWLF/nxxx9rdN/AwED69OlD586dcXJy4tChQ6xfv55XXnmlTm0SBEEwGFlW1rv49S0oyQEjU2XeRa9/gLGp7rKM4gx2JO0gKimKQ2mHqgWLQOdA1L5qwhyeoHLau5RcUAKGz5o1WLS5MXk/r6SciWsOcjAxG0tTY5b/rTMhAXeeuN/QRMhoht599907Hl+9ejV9+vSp0z3d3d35448/eOutt1i0aBElJSUEBwezZcsWhgwZUqN7TJ06lc2bN6PRaCgtLcXX15cPP/yQ6dOn16lNgiAIBpF3FX5+Dc5uU8oeHeEv/wb39gBcK7rGjuQdaBI1HE47jIyse2mQcxBqlZoBvgPwtvWmMieHpAkTKI07jbGTEz5rVmNx09OBmQWlvLD6ACcv52FrYcKa8Y/T2depQd/uvUiyLN//qoeAJEl2QG5ubi52dnZ3vKaoqIjTp0/Trl07rKya1vrwQu1UfdaJiYkkJCQwePBgOnbsaOhmCYLQlMkyHPsWts+EklwwNoPQt6DnNNJLs4lKikKTqOFo+tFqwaKDSwelx0IVhpeNl+54RXY2yRMmUnr6NMbOzviuWY35TY/3p+aW8Pyq/ZxLL8DZ2oz/TOxKe097vby1vLw87O3tAexlWc6r6etET4YgCIIgPKjcy8p27Oeub4Pg+Rhpgz4kquAiUZpJtwWLYNdgJVj4huFp43nb7Sqys0keP4HSM2eUgBG5BvObNrFMziziuVV/kpJVjIe9Bd9M6oa/q43e32ZtiZAhCIIgCHUly3D0v7D9bSjNI9XMkqjgoWiMSji2a0q1Szu5dkKtUoJFC+u7T8qsyM4medx4SuPjMXZxUXowbgoYZ9PyeX7lftLzS1E5W/HNpG60dGycve8iZAiCIAhCXeSkwJapXE38nShrKzQt/DluVA6Z+3WXPOr2KGpfZY7FvYJFlYqsLCVgnD2LsasLvpGRmLdqpTsfeymHv319gJyictq42/LfSV1xszX8olt3I0KGIAiCINSGLHNl35dEHfoSjbkxsT5V8yjKkZCUYKFSM8BnAO7W7jW+bUVmphIwEhLuGDD+OJfBi/89TEFpBR29HYgc/zgOVndeZ6OxECFDEARBEGrgUv4los78D03cWk5SBvbKEIWExGPuj+l6LNys3Gp9byVgjKM04Rwmrq74REZi3krZfqGgtIIF287wnz+TkGXo3sqZFS90wca88f8Kb/wtFARBEAQDSclP0T0VcirzlO64kSzT2coLdYcXGKBS42LpUuc6KjIySBo3jrJz5zFxc8Mncg3m1/d3+u10GrM3neRqrrJL9jNdWjJ3eBAWDbC9RH0QIUMQBEEQbpKSl8L2pO1oEjWczjqtO24kyzxeUkqYeQv6D1yIi2fNdky9l4pr10gaN56y8+cxcXfHN3INZioV1/JLeX/LKX6OvQqAj5MVEX/tQM/WdQ8zhiBChiAIgtDsJeUloUnUoEnScCbrxorIRkg8XlKKuqCA/mUyzv3eg8cnQx1XR75ZxbVrJL0wjrILFzBp0QLfyDWY+vjwv0MpfLT1NLnF5RgbSUwK8eO1/o9gadY0ei9uJkKGIAiC0CxdzL2oGwqJz76xwaOxZExX5yDU6cn0SzmBk1YLqhAY9iU4tbrHHWuuPD2d5HHjlYDh4YFv5BquWjszc+V+/jifCUCQlx0f/zWYIC/9LLDVEETIEARBEJqNC7kXdD0WCdkJuuMmkgndPLoR5tOffteScYxeABXFYGoNg96HLhPrpfcCoDwtneQXXqAsMRETDw9arl7N1xfL+VfUbkortFiYGvF62CNM6OmHiXH91GkoTSZkSJLkAOxAabMJ8IUsyysM2ypBEAShsTufc14XLM7lnNMdN5FM6ObZjYG+A+nr3ReHggz46WVI+VO5wK83DFsEjqp6a0u1gOHpQeknixm5OZlTV5SVunu2dmbeiA74OlvXW52G1GRCBpAP9JZluUiSJGvgpCRJG2RZzjR0w4Qbxo0bR3R0NImJiYZuiiAIzZQsy5zLOacbCjmfe153zsTIhO4e3QnzDaOfTz/sze1BWwl/LoGdH0JFCZjZgPpD6DwOJKne2lWelkby316gLCkJY09Pto5/hy83XEQrg72lKbOHtGNk55ZI9VinoTWZkCHLciVQdL1oDkjX/wm1cOrUKSIiIti1axcZGRk4OzvTt29fZs2aRfv27Q3dvFrZsGEDP/zwAwcPHiQ1NRVvb2+GDh3KO++8g4ODg6GbJwhCA5JlmYScBF2PxcXci7pzJkYm9PDsgdpXTR/vPkqwqHLtLPz0Elw6qJRb9VXmXjj41Gv7ylNTSXrhBcqTkql0a8Gsnn/n6KlCAJ7q6Mm7QwNxtTWv1zobg3oLGZIk9QamA50BD2CELMubbrnm5evXtACOA6/KsnygFnU4AL8DAcB0WZYz6qn5zcKGDRsYPXo0Tk5OTJw4ET8/PxITE1m1ahXr16/n+++/Z8SIEYZuZo29+OKLeHp68vzzz+Pj48OJEydYvHgxv/zyC0eOHMHS0tLQTRQEQY9kWeZs9lk0SRo0iRoS8xJ150yNTOnp2RO1Sk2odyh2Zrfsvl1ZAfsWw655UFkK5nZK78Vjf6vX3guA8qtXSXphHOXJyeQ7uPJKpwmkl1viaW/BhyOC6Ne25quCNjX12ZNhjRIcvgY23HpSkqRngc+BKcB+4DVguyRJbWRZTr9+zbG7tEkty/IVWZZzgI6SJLkDGyRJWi/Lclo9voeH1vnz5xk7diytWrVi9+7duLq66s5NmzaNkJAQxo4dS2xsLK1a3X32dGFhIdbWDTNWKMsyJSUldw0L69evp0+fPtWOde7cmRdeeIG1a9cyadKkBmilIAgNSZZl4rPjdT0WSXlJunNmRmb08LrRY2FrZnvnm6Sfhk0vwZUjSrn1AHjqC7BvWe/tLb9yRQkYKSmk2zjzRtcXybB2ZFx3FW8MbNMkVu18EPU2bVWW5V9lWZ4ty/LGu1zyOrBCluXVsizHoYSNImDCTffoJMty0B3+XbmlrjSUQBNyt/ZIkmQuSZJd1T/gLv+3NQ8LFiygqKiI5cuXVwsYAC4uLixbtozCwkLmz5+vOz5nzhwkSSIuLo4xY8bg6OhIr169dOc3bdpEUFAQFhYWBAUFsXHjnT96rVbLwoULad++PRYWFri7uxMeHk52dna161QqFUOHDmX79u106dIFS0tLli1bdtf3dGvAAHQ9MadPn77tnCAITZMsy8RlxrHw8EKGbhzKqC2jWHFiBUl5SZgZmdHPux8fh3zM78/+zqJ+i3jK/6k7B4zKCtjzGSzrrQQMc3sYvgSeW6+3gHH++b9RnpLCVStn3ugxBUc/b378ew/mDGv/0AcMaKA5GZIkmaEMo0RUHZNlWStJ0g6gew3v4Q4UybKcL0mSPdAb+Pc9XjITeK/urX64bNmyBZVKRUjInXNZ7969UalUbN269bZzo0aNIiAggHnz5iHLMgAajYann36awMBAIiIiyMzMZPz48bRsefs3anh4OGvWrGH8+PFMnTqVixcvsnjxYo4ePcrevXsxNTXVXRsfH8/o0aMJDw9n8uTJtGnTplbvMzU1FVCCkyAITVdVsKgaCrlUcEl3ztzYnBCvENQqNb1b9sbatAa9q2mnlN6Lq8eUcsBAeGoh2Hnqpf0lKZc4PXosFhmpXLF25p3eLzF26ONMCfXHzKRpP5ZaGw0Vo1wAY+DWoY00oG0N7+ELLJeUabcSsEiW5RP3uD4CZXimii1w6S7X3p0sQ3nR/a9rSKZWtRozzM3N5cqVKwwfPvye1wUHB7N582by8/Oxtb3xV0DHjh359ttvq1371ltv4e7uTkxMDPb2yiSq0NBQ1Go1vr6+uutiYmJYuXIla9euZcyYMbrjffv2ZdCgQaxbt67a8XPnzrFt2zYGDhxY4/d3s08++QRjY2NGjhxZp9cLgmA4sixzKvOUbijkcsFl3TkLYwtCWoag9lWChZWpVc1uWlkOMf+C3+eDthws7GHwfAh+tt7nXlQ5czSezPBJOOVlcNnahW+ffYv/jAultVvz61BvMn011yeIdqrF9aVAaVW5zo8ElRfBPP0k3TqbdQXMaj4vIj8/H6BacLiTqvN5eXnVrp0yZUq1665evcqxY8eYMWOGLmAAhIWFERgYSGFhoe7YunXrsLe3JywsjIyMG/N0O3fujI2NDbt27aoWMvz8/OocML799ltWrVrFm2++SUBAQJ3uIQhCw5JlmRMZJ9AkaohKiuJK4Y3RcUsTS12PRYhXSM2DRZXUE0rvRWqsUm7zJAz9F9i2qMd3cENJeSVfr99L4IIZuBdlc8XGlcwP/sXKgY9hZNQ8H4ZsqJCRAVQCt06hdQdSG6gNzVZVYKgKG3dztzDid303wCpJScpEqzv9Im/Tpg1HjhzRlRMSEsjNzcXN7c5bH6enp9+zrpras2cPEydOZODAgXz00Ud1uocgCA1DK2uJvRZLVFIUUUlRXC28qjtnaWJJaMtQwnzD6OXVq/bBAqCiTJl7sedT0FaApSMMXgAdRuqt9+LAxSwWrNnJy5s/w704h2ynFrSLjMQzoH4fhW1qGiRkyLJcJknSYaA/sAlAkiSj6+XFDdGGOjO1UnoOGpNaftPZ29vj4eFBbGzsPa+LjY3Fy8sLO7vqj3o9yKOgWq0WNzc31q5de8fzt05CrUtdx48fZ9iwYQQFBbF+/XpMTJpMB50gNBtVwWJ74naikqJIK7oxem5pYkmfln1Qq9T09OqJpckDPH5+9bjSe5F2Uim3HQpDPgdb/TwmmldSzse/nmHnzqN8EvNv3IpzKPdoSZfvv8HM/eF9NLWm6nOdDBug9U2H/CRJ6gRkybKcjDI/IlKSpEPAAZRHWK2B1fXVBr2QpFoNTTRWQ4cOZcWKFcTExFR7QqTKnj17SExMJDw8/L73qppzkZCQcNu5+Pj4amV/f3927NhBz5499bJuxfnz5xk0aBBubm788ssv2NjY1HsdgiDUjVbWciz9GJokZSgkvehGz6WViRV9vK8HC8+eWJhYPFhlFaWwewHs+RzkSrByhicXQPu/6q33YvupVN7ZdBLj1MvMj/k3rsW5mKhUtP5PJKZ36b1tburzT74uwK6bylWTLiOBcbIs/yBJkiswF2UxrmPAILHORcOYPn0633zzDeHh4ezevRtnZ2fduaysLKZMmYKVlRXTp0+/7708PDzo1KkTkZGR1eZlREVFERcXV23i5zPPPMOSJUv44IMPmDdvXrX7VFRUUFBQUOfVOVNTU1Gr1RgZGbF9+/bbekUEQWh4ldpKjl07hiZRw46kHaQX3wgWNqY29PHuQ5hvGD29emJuXE8rXF4+ouw5kh6nlAP/Ak9+Cjb6+ZmQUVDKez+dYuuJq3gUZPDZvmU4Fudi1qoVvpFrMBE/i3TqLWTIshzNfZb5lmV5MY19eOQhFRAQQGRkJM899xwdOnS4bcXPjIwMvvvuO/z9/Wt0v4iICIYMGUKvXr2YMGECWVlZLFq0iPbt21NQUKC7LjQ0lPDwcCIiIjh27BhqtRpTU1MSEhJYt24dX3zxRZ2fBBk0aBAXLlzgzTffJCYmhpiYGN05d3d3wsLC6nRfQRBqp1JbyZH0I0qwSN5BRvGNSd42pjb09e6LWqWmh2cPzIzN6q/iilKI/hj2fnG998IFhnwK7fWzcrEsy/wce5X3Np8iq7AM76IMvji4AsvCbMz8/fFds1oEjFuIwetmZNSoUbRt25aIiAhdsLh575KgoKAa36vq8dPZs2czc+ZM/P39Wb16NT/99BPR0dHVrl26dCmdO3dm2bJlzJo1CxMTE1QqFc8//zw9e/as8/s5fvw4QLUFxKqEhoaKkCEIelSpreRw2mE0SUqPRWbJjb0qbU1t6evTF7Wvmu6e3es3WFS5dFjZc+TaGaUc9LQyudPa+d6vq6P0/BLe2XSS7aeUzvdeFkXM2rMKKTcTs9b++K5Zg4lYn+c2UtXiSg+766t+5ubm5t42sbFKUVERp0+fpl27dlhZ1WFGs9BkVH3WiYmJJCQkMHjwYDp27GjoZglCo1ahrVCCxfUei6ySLN05WzNb+nn3Q61S092jO6bGpve40wMoL4HoefDHIpC1YO0GQz+Hdk/ppTpZltl07DJzNseRW1yOiZHEm4EW9Fk2h8pr1zAPaI3PmjWYOOsn3DQWeXl5VUPj9rIs59X0daInQxAEQbirCm0FB1MPoknSsDN5Z7VgYWdmR3+f/qhVarq16Ka/YFEl5YAy9yLjrFIOfhYGfQxWTnqpLjW3hLc3nuC3M8q8kvaedszvao/Z9FeouHYN84AAfNasfugDxoMQIUMQBEGoplxbrgSLRCVYZJfe2GfI3txeCRa+arp6dMXUSM/BAqC8GHZ+CPu+AmSwcYehC6Htk3qpTpZl1h2+xAc/x5FfUoGpscS0/gGM95a4PGE8FdcyMH/kESVgOOkn4DwsRMgQBEEQKNeWc+DqAV2PRU5pju6cg7mDrsfi8RaPN0ywqJK0T+m9yDqvlDuOhoHz9NZ7cSWnmJkbTvD72WsABLe0Z8HIjqgK00l6YRyVGRmYt2mDz+qvRcCoAREyBEEQmqnyynL+vPqnLljkld0Yanc0d6S/r9Jj8XiLxzExauBfF2VFsPMD+PPfgAy2Hsp27I/UbduB+5Flme8PpvDR1tMUlFZgZmLEPwY8wuQQPyovXiBp3HglYLRtqwQMR0e9tONhI0KGIAhCM1JeWc6+q/uUoZCUneSX3dhuwMnCiQE+A1Cr1HR279zwwaJK4l6l9yL7olJ+9HlQfwSWdVtT535SsoqYueEEMeeUR28f9XFgwciOtHazoTQhQQkYmZmYt2uHz9erRMCoBREyBEEQHnJllWXsu7IPTZKGXcm7yC+/ESycLZwZ4DsAta8SLIyNjA3X0NIC+O19OLBcKdt5wVNfQsAAvVSn1cp8sz+Jj389Q1FZJeYmRkwf2IbxPf0wNpIoOXuW5HHjqczKwjywHT6rRMCoLREyBEEQHkKllaX8cfkPopKi2JWyi4LyG4vkuVi66HosHnN7zLDBosrF3fDTK5CjbMDIYy+A+gNla3Y9SMos5M31sey/qDwt01XlxCcjg/FzUbaRKIk/S/J4JWBYBAbi8/UqjOu4OnFzJkKGIAjCQ6K0spS9l/eiSdIQnRJNYXmh7pyrpauux+JRt0cbR7AAKM2HqPfg0CqlbO8Nw74E/356qU6rlVnzRyILtsdTXF6Jpakxbw1qw9+6q3TbsZfExys9GNnZWLRvj8+qlSJg1JEIGYIgCE1YSUWJLlj8fun3asHCzdKNMFUYal81ndw6YSQZGbCld3AhGn56FXKTlXKXCRA2F8xt9VPdtQLeXB/LoSTlkdzurZz55OlgfJxvLL5YcuaMEjBycrAIClIChr1+elOaAxEyBEEQmpiSihJiLsegSVSCRVFFke6cu5U7Yb5hqFVqOrp2bHzBAqAkD6LegcNrlLKDDwxbDK1C9VJdpVZmVcwFPtOcpbRCi7WZMTOfbMeYrj663guAktOnSR4/QQkYHTooAeMuK0QLNSNChiAIQhNQXFHMnkt7iEqK4vdLv1NcUaw718K6hRIsfNUEuwY3zmBR5dxvsHkq5F1Syl1fhP7vgbmNfqpLz+eNdbEcS1HW/QgJcCHirx1o6Vh964iSuDglYOTmYhEcjM/KFSJg1AMRMgRBEBqpovIi9lzegyZRw57Le6oFC09rT12PRQeXDkjSPTfBNrySXNj+Nhz9r1J2VMHwr0DVSy/VVVRqWb7nAgt3JFBWocXW3IS3h7Tj2ce9b/taFZ86RfKEiWhzc7HoGIzPypUY2+pnyKa5ESFDqFfjxo0jOjqaxMREQzdFEJqkovIidl/ajSanvLlJAAAgAElEQVRJw55LeyipLNGd87LxQu2rJsw3jCCXoMYfLKqc1cCWaZB/BZCg2xTo/w6YWeuluvjUfKavP07spVwA+rRxJeKvHfCwt7zt2psDhmXHjnivXCECRj0SIaOZOXXqFBEREezateu2rd7bt29v6ObVSnx8PEuXLmX//v0cOXKE0tJSLl68iEqlMnTTBKFWisqL+P3S72gSNcRcjrk9WKjUDPQdSKBzYNMJFgDF2UrvxbG1StnJX+m98O2ul+rKK7X8O/o8i3YmUF4pY2dhwrtPtefpx7zu+HUrPnmK5AkT0OblYdmpkxIwbPQzbNNciZDRjGzYsIHRo0fj5OTExIkT8fPzIzExkVWrVrF+/Xq+//57RowYYehm1ti+ffv48ssvCQwMpF27dhw7dszQTRKEGissLyQ6JZqopChiLsdQWlmqO+dt6630WKjCCHRqYsGiSvyvsOU1KEgFJOj+MvR9G8ys7vvSujh1JZfp62KJu6osjT6gnTsfjQjC3c7ijtcXnzhB8sRJSsB49FG8V6zA2EY/PSvNmQgZzcT58+cZO3YsrVq1Yvfu3bi6uurOTZs2jZCQEMaOHUtsbCytWrW6630KCwuxtm6Yb0RZlikpKcHS8vYuToBhw4aRk5ODra0tn376qQgZQqNXUFZA9KVoNIka9l7eS5m2THfOx9YHtUqN2ldNW6e2TTNYABRlwbaZEPu9UnYOUHovfLrppbqyCi2LdyawJPo8FVoZBytT3h/WnmEdPe/6NSyOjVUCRn4+lo89hvfy5SJg6EkjnoIs1KcFCxZQVFTE8uXLqwUMABcXF5YtW0ZhYSHz58/XHZ8zZw6SJBEXF8eYMWNwdHSkV68bk7Q2bdpEUFAQFhYWBAUFsXHjxjvWrdVqWbhwIe3bt8fCwgJ3d3fCw8PJzs6udp1KpWLo0KFs376dLl26YGlpybJly+76npycnLAVY6dCI5dfls+W81t4deerhP4Qysw9M9mVsosybRkqOxWTO0xm/VPr+XnEz0x7bBrtnNs13YBx+mdY8oQSMCQj6DEVpuzRW8A4cSmXYYtj+HLnOSq0MoODWhD1j1CGd7rz8AhA8fHjyhyM/Hwsu3QWAUPPRE9GM7FlyxZUKhUhISF3PN+7d29UKhVbt2697dyoUaMICAhg3rx5yLIMgEaj4emnnyYwMJCIiAgyMzMZP348LVu2vO314eHhrFmzhvHjxzN16lQuXrzI4sWLOXr0KHv37sXU9Ma20fHx8YwePZrw8HAmT55MmzZt6ukrIAgNJ68sj+gUpcfijyt/UK4t153zs/fTTd58xPGRphsoblaYCb++CSfXK2WXNvCXJdCyi16qKymv5MvfEli2+wKVWhlnazPmDg9iSLDHPV9XfOwYyZMmoy0owKpLF7yXLcWogXpmmysRMu5DluVqj401BpYmlrX6wZSbm8uVK1cYPnz4Pa8LDg5m8+bN5OfnV+sh6NixI99++221a9966y3c3d2JiYnB/vpqeKGhoajVanx9fXXXxcTEsHLlStauXcuYMWN0x/v27cugQYNYt25dtePnzp1j27ZtDByon+2cBUFfcktz2ZWyC02ihn1X91GhrdCda2XfSjcU0tqh9cMRLKrEbYatr0PhNaX3oudrEPoWmN55LsSDOpqczfT1sZxLV/ZieaqjJ3OeCsTZxvyerys6epSUSZPRFhZi9fjjeC/9twgYDUCEjPsoriim27f66eqrq/1j9mNlWvPJU/n5yo6L9xtaqDqfl5dX7dopU6ZUu+7q1ascO3aMGTNm6AIGQFhYGIGBgRQW3ljWeN26ddjb2xMWFkZGRobueOfOnbGxsWHXrl3VQoafn58IGEKTkVuay87knWiSNPx59c9qwaK1Q2vUvmrUKjX+Dv4GbKWeFGbAL2/AqevDpK7t4C9fgVdnvVRXUl7J51FnWbnnAloZXGzM+fAvQQwKanHf1xYdOUrK5OsBo2tXJWBY6WcCqlCdCBnNQFVgqAobd3O3MOLn51etnJSk7JIYEBBw2z3atGnDkSNHdOWEhARyc3Nxc3O7Y53p6en3rEsQGpuckhx2puxEk6hh/9X9VMg3gkWAY4ASLHzVtHK4+wTqJu/URtj6TyjKBMkYQl6H3tPB5N69CXV1KDGLN9fHciFD+QNmxKNevPdUIA5WZvd9bdGRI0oPRlERVt264f3vJSJgNCARMu7D0sSS/WP2G7oZ1Via3Plpi7uxt7fHw8OD2NjYe14XGxuLl5cXdrcspXu3pztqQqvV4ubmxtq1a+94/tZJqA9SlyDoS3ZJNr8l/4YmUcOB1ANUypW6c484PqLrsfCzf8hDckG6Ei5Ob1bK7kHKkyOenfRSXVFZBQu2x7Pmj0RkGdztzJk3ogP927nX7PWHD5My+UUlYDzxhBIwxM+YBiVCxn1IklSroYnGaujQoaxYsYKYmJhqT4hU2bNnD4mJiYSHh9/3XlVzLhISEm47Fx8fX63s7+/Pjh076NmzpwgQQpOSVZKlCxYHUw9WCxZtndrqJm+q7FWGa2RDkWU4+SP8Mh2Ks8DIBELegJB/gsn9exPq4s8Lmbz1YyxJmcrmb6M6t2T20EDsLU3v80pF0aFDJL8YjlxUhHWP7rT86isRMAxAhIxmYvr06XzzzTeEh4eze/dunJ2ddeeysrKYMmUKVlZWTJ8+/b738vDwoFOnTkRGRlablxEVFUVcXFy1iZ/PPPMMS5Ys4YMPPmDevHnV7lNRUUFBQQEODg719C4F4cFkFmfeCBZpB9HKWt25dk7tdJM3fex8DNjKBpafpkzsPPOzUm7RAYYvAY9gvVRXWFrBJ9vO8J99yrCsh70FEX/tQJ82dx5yvZOigwdJDp9yPWD0oOWSrzCy0M9EVOHeRMhoJgICAoiMjOS5556jQ4cOt634mZGRwXfffYe/f80mqEVERDBkyBB69erFhAkTyMrKYtGiRbRv356CggLddaGhoYSHhxMREcGxY8dQq9WYmpqSkJDAunXr+OKLLxg5cmSd3lNubi6LFi0CYO/evQAsXrwYBwcHHBwceOWVV+p0X6F5ySjOYEfSDqKSojiUdqhasAh0DtTNsfC28zZgKw1AliH2f8qjqSU5YGQKoW9Cr3+Acc16E2pr77kM3voxlkvZyhN9o7v6MOvJttha1Ly+wgMHSAmfglxcjHXPnrT8arEIGAYkQkYzMmrUKNq2bUtERIQuWNy8d0lQUFCN71X1+Ons2bOZOXMm/v7+rF69mp9++ono6Ohq1y5dupTOnTuzbNkyZs2ahYmJCSqViueff56ePXvW+f1kZ2fzzjvvVDv22WefAcqQjggZwt1cK7rGjuQdaBI1HE47jIysOxfkHIRapQyFtLS9fd2XZiHvKvz8GpzdppQ9Oiq9Fy1q/jOiNvJLypn3yxm+O5AMgJeDJfNHBtOztUut7lO4/wApU64HjF69aLl4kQgYBiZVLa70sJMkyQ7Izc3NvW1iY5WioiJOnz5Nu3btsBKzjx9qVZ91YmIiCQkJDB48mI4dOxq6WYIepRelE5UUhSZRw9H0o9WCRQeXDrq9QrxsvAzYSgOTZTj+HWyboWzNbmymrHnRc5reei9+P3uNmT/GciVX2RTub919eXNQW2zMa/c3cOGff5Iy5e/IJSVYh4QoAcNcP0+7NEd5eXlVQ+P2sizn1fR1TaYnQ5KkNsAPNx1qA4yWZXmTgZokCEIjl1aYRlRSFFFJUbcFi2DXYN3kTU8bTwO2spHIvaxsx34uSil7Pqas2unWTj/VFZfz4c9xrDt8CQBfZys+eTqYJ1o53+eVtyvct4+Uv7+kBIzeIbRcJAJGY9FkQoYsy/FAJwBJkmyARCDKkG0SBKHxSS1M1fVYHLtWfdO8Tq6ddEMhLazvv4hTsyDLcPS/ypbspXlK70XfWdD9VTDWz6+I306nMWvjCdLySpEkGN/DjzcGPoKVWe3rK/zjDyVglJZiExqK16IvMTLTzxMvQu01mZBxi2HAb7IsF973SkEQHnpXC66iSdKgSdIQe636ejCPuj2K2lfNAN8BIljcKicFtkyF8zuVslcXpffCVT97BuUUlTF3Sxwbjl4GoJWLNfNHBtNF5VSn+xXs3cull15WAkafPnh9+YUIGI1MvYUMSZJ6A9OBzoAHMOLWoQxJkl6+fk0L4DjwqizLB+pQ3TPAfx6sxYIgNGWXCy6zI0mZvBmbcSNYSEhKsFCpGeAzAHfrmi3c1KzIMhxeA5p3oCwfTCyg79vQ/WUwMtZLldtPpTJ700mu5ZdiJMGkkFa8HvYIFqZ1q68gZi+XXnoJuawMm7598fpioQgYjVB99mRYowSHr4ENt56UJOlZ4HNgCrAfeA3YLklSG1mW069fc+wubVLLsnzl+jV2QA/g/+qx7YIgNAGX8i/phkJOZp7UHZeQ6OzemTDfMAb4DsDNquZrKjQ72UlK78WFaKXs3U1ZtdPl9m0C6kNWYRnvbT7FluNXAGjtZsOCkcE86uNY53sW7Inh0ssvKwGjf39a/utzJBEwGqV6CxmyLP8K/ArcbYfB14EVsiyvvn7NFGAIMAH4+Po9arI27XBAI8tyyb0ukiTJHLh55s+9dwcTBKFRSslPQZOoDIXEZcbpjhtJRnR276wbCnGxrN3jjs2OVguHv4ao96CsAEwsof870G2K3novfjlxlXc2nSSzsAxjI4nw3q2Y2j+gzr0XAAW7d3PplVeVgDGgPy0/FwGjMWuQORmSJJmhDKNEVB2TZVkrSdIOoHstb/cMsLwG180E3qvlvQVBaARS8lLYnrQdTaKG01mndceNJCMed3+cMN8w+vv2F8GiprIuwuZXIXGPUvbprvReOOtnd9iMglLe++kUW09cBaCNuy2fjupIh5b293nlvRX8/rsSMMrLsQ0bgNdnn4mA0cg11MRPF8AYSLvleBrQtqY3kSTJHugKPF2DyyNQhmeq2AKXalqXIAgNKykvSddjcSbrjO64kWTE4y0eR+2rpr9Pf5wta/+IY7Ol1cLBlbDjPSgvAlMr6P8edH0RjIzqvTpZlvk59irvbT5FVmEZJkYSL/Xx55V+AZiZPFh9+dHRXH516vWAEYbX558hmepn7Q6h/jSpp0tkWc4FajSLS5blUqC0qnyXIRxBEAzoYu5F3RyL+Owbm+sZS8Z0bdEVtUpNP59+OFnU7emDZi3zvNJ7kaQsuY9vLxi+CJz0swV9en4J72w6yfZTyt+S7TzsWDAymCCvB+u9AMjfuYtL06ZBeTm2ajVen30qAkYT0VAhIwOo5PaA4A6kNlAbBEFoBC7kXtD1WCRk39jJ10QyoZtHN8J8w+jn0w9Hi7pPDGzWtJWwfxn8NhcqisHUGsLehy4T9dZ7sfn4Fd7bfIqconJMjCRe7RfA3/v4P3DvBUD+zp1cmvaaEjAGDcJrwXwRMJqQBgkZsiyXSZJ0GOgPbAKQJMnoenlxQ7RBEATDOZ9zXhcszuWc0x03kUzo5tmNgb4D6evdFwcLsSPvA8k4Bz+9DCl/KmVVCAxfDI4qvVSXnlfCrI0n2XFa6b1o72nHgpEdCfS889YNtZX/229ceu0fSsAYPAivBQuQTJpUB3yzV5/rZNgArW865CdJUicgS5blZJT5EZGSJB0CDqA8wmoNrK6vNgiGN27cOKKjo0lMTDR0UwQDkmWZcznndEMh53PP686ZGJnQ3aO7rsfC3vzBu9ObPW0l/LkEdn4IFSVgZgNhc6HzeL31Xmw8epk5m0+RV1KBqbHEtP4BhIf6Y2pcP/XlRUVx+R+vQ0UFdk8+ief8T0TAaILq8xPrAuy6qVw16TISGCfL8g+SJLkCc1EW4zoGDJJl+dbJoIIenTp1ioiICHbt2nXbLqzt27c3dPNqZePGjSxdupQTJ06QmZmJq6srTzzxBHPmzKnVjrJC/ZBlmYScBF2PxcXci7pzJkYm9PDsgdpXTR/vPiJY1KdrZ+Gnl+DSQaXcqg8MWwQOPnqpLjW3hFkbT7DzTDoAHbzs+XRUR9q0qL9VAvI0Gi6//k8lYAwZgucnH4uA0UTV5zoZ0cA9Z1fKsrwYMTxiMBs2bGD06NE4OTkxceJE/Pz8SExMZNWqVaxfv57vv/+eESNGGLqZNXbixAkcHR2ZNm0aLi4upKam8vXXX9O1a1f27dsndlVtALIsczb7LNsTtxOVFEViXqLunKmRKT09e6JWqQn1DsXOrH660IXrKitg3yLYFQGVpWBmCwM/gsf+BnqY6C7LMusPX2Luz3Hkl1RgZmzEa2EBvBjSCpN66r0AyNuu4fI/rweMoUPx/DhCBIwmTHxyzcT58+cZO3YsrVq1Yvfu3bi6uurOTZs2jZCQEMaOHUtsbCytWt199nlhYSHW1tYN0WRkWaakpARLS8s7nn/33XdvOzZp0iRatmzJv//9b5YuXarvJjZLsixzJuuMMhSSpCEpL0l3zszIjB5eN3osbM3EGnh6kX4aNr0EV44o5dYD4KkvwL6lXqq7klPMrI0niI6/BkBHbwc+HRlMgHv9fr5527Zx+Z9vQGUldsOewjMiAslYPwuFCQ2j/gfrhEZpwYIFFBUVsXz58moBA8DFxYVly5ZRWFjI/PnzdcfnzJmDJEnExcUxZswYHB0d6dWrl+78pk2bCAoKwsLCgqCgIDZu3HjHurVaLQsXLqR9+/ZYWFjg7u5OeHg42dnZ1a5TqVQMHTqU7du306VLFywtLVm2bFmt3qebmxtWVlbk5OTU6nXCvcmyTFxmHAsPL2ToxqE88/MzrDixgqS8JMyMzOjn3Y+PQz7m92d/Z1G/RTzl/5QIGPpQWQF7PoNlvZWAYW6vLKr13Hq9BAxZlvnhYDID/7Wb6PhrmJkYMXNwW36c0r3+A8avv+oChv3wYSJgPCRET0YzsWXLFlQqFSEhIXc837t3b1QqFVu3br3t3KhRowgICGDevHnIsgyARqPh6aefJjAwkIiICDIzMxk/fjwtW97+gy48PJw1a9Ywfvx4pk6dysWLF1m8eDFHjx5l7969mN70OFp8fDyjR48mPDycyZMn06bN/XeDzMnJoby8nNTUVBYuXEheXh79+/ev6ZdGuIuqYLE9aTtRiVFcKrixlp25sTm9vHqh9lWGQqxNG6Z3q1lLO6X0Xly9vn19wEB4aiHYeeqluss5xcz4MZY9CRkAPOrjwIKRHWntZlPvdeX98guXp795PWAMx2PeRyJgPCREyLgPWZaRi4sN3YxqJEvLWi0ulpuby5UrVxg+fPg9rwsODmbz5s3k5+dja3vjr5SOHTvy7bffVrv2rbfewt3dnZiYGOztlUl8oaGhqNVqfH19ddfFxMSwcuVK1q5dy5gxY3TH+/bty6BBg1i3bl214+fOnWPbtm0MHDiwxu/viSeeID5eWcjJxsaG2bNnM3HixBq/XrhBlmVOZpzUDYVcLrisO2dhbEFIyxDUvmp6t+yNlamVAVvajFSWQ8y/4Pf5oC0HC3sY9Al0/D+9zb349kAyEb+coaC0AnMTI6YPbMP4nn4YG9V/fblbt3Jl+pug1WI/YgQeH34gAsZDRISM+5CLi4l/rLOhm1FNmyOHkaxq/gM+Pz8foFpwuJOq83l5edWunTJlSrXrrl69yrFjx5gxY4YuYACEhYURGBhIYWGh7ti6deuwt7cnLCyMjIwM3fHOnTtjY2PDrl27qoUMPz+/WgUMgNWrV5OXl8eFCxdYvXo1xcXFVFZWYqSHR/ceRrIscyLjBJpEDVFJUVwpvKI7pwsWKjW9vUSwaHCpJ5Tei9TrW9m3eRKG/gtsW+ilupSsImZsiGXvuUwAuvg6Mn9kMK1c67/3AiB3y89ceestJWD89a94fDBXBIyHjAgZzUBVYKgKG3dztzDi5+dXrZyUpEz0Cwi4fWvoNm3acOTIEV05ISGB3Nxc3NzuvPV2enr6Peuqie7db+yx93//93+0a9cOgE8//bTW92outLKW2GuxRCVFEZUUxdXCq7pzliaW9G7ZG7Wvml5evUSwMISKMmXuxZ5PQVsBlo4weAF0GKmX3gutVmbtgWQifjlNUVklFqZGvDmwLS/0UOml9wIgd8sWrrw1QwkYI5/GY+5cJPGHwUNHhIz7kCwtaXPksKGbUY10l6ct7sbe3h4PDw9iY2PveV1sbCxeXl7Y2VV/1PBuT3fUhFarxc3NjbVr197x/K2TUB+kLgBHR0f69evH2rVrRci4RVWwqHrcNK3oxhI1liaWhLYMRa1SgoWlyYN9DsIDuHJMWbUz7aRSbjsUhnwOtjXatqnWkjOLePPH4/x5IQuArion5o8MRuWiv3k2uZs3c2XGTNBqcRg1khbvvy8CxkNKhIz7kCSpVkMTjdXQoUNZsWIFMTEx1Z4QqbJnzx4SExMJDw+/772q5lwkJCTcdq5qbkQVf39/duzYQc+ePR84QNRUcXExubm5DVJXY6eVtRxLP4YmSRkKSS+60XNkZWJFqHcoA30H0tOrJxYmFgZsqUBFKexeAHs+B7kSLJ1gyKfQ/q966734z75EPtkWT3F5JZamxswY3JaxT/hipKfeC4CcTZu4OnMWyDIOo0bR4v05ImA8xETIaCamT5/ON998Q3h4OLt378bZ+cZ22VlZWUyZMgUrKyumT59+33t5eHjQqVMnIiMjq83LiIqKIi4urtrEz2eeeYYlS5bwwQcfMG/evGr3qaiooKCgAAeHuu1XkZ6eftswTGJiIr/99htdunSp0z0fBpXaSo5dO4YmUcOOpB2kF98IFtam1vTx7kOYbxg9PUWwaDQuH1HmXlw7rZQD/wJPfgo2rvd+XR0lZhTy5o+xHLio9F480cqJ+U93xMdZv39Q5WzcxNVZ1wPGs8/S4r13RcB4yImQ0UwEBAQQGRnJc889R4cOHW5b8TMjI4PvvvsOf3//Gt0vIiKCIUOG0KtXLyZMmEBWVhaLFi2iffv2FBQU6K4LDQ0lPDyciIgIjh07hlqtxtTUlISEBNatW8cXX3zByJEj6/SeOnToQP/+/enUqROOjo4kJCSwatUqysvL+fjjj+t0z6aqUlvJkfQjSrBI3kFG8Y1JtjamNvTx7oPaV00Prx6YG5sbsKVCNeUl8PvHsPdLpffCyuV674V+Vt7VamXW/JHI/O1nKCnXYmVmzMwn2/FcVx+99l4A5Py4gauzZysB4/+epcW7ImA0ByJkNCOjRo2ibdu2RERE6ILFzXuX1Ga/j6rHT2fPns3MmTPx9/dn9erV/PTTT0RHR1e7dunSpXTu3Jlly5Yxa9YsTExMUKlUPP/88/Ts2bPO7+fvf/87W7duZdu2beTn5+Pm5oZarWbWrFl06NChzvdtKiq1lRxOO4wmSemxyCzJ1J2zNbWlr09f1L5qunt2x8zYzIAtFe7o0iGl9yLj+hBj0NMweD5Yu+ilugvXCnhzfSyHkpRF8Hr4O/PJ08F4O+l/ODjnxx+5OvsdkGUcx4zG/Z13avUYvtB0SVWLKz3sJEmyA3Jzc3Nvm9hYpaioiNOnT9OuXTusHoJ5GMLdVX3WiYmJJCQkMHjw4Cax10mFtkIJFtd7LLJKsnTnbM1s6evdl4GqgTzh8YQIFo1VeTHsmgf7FoOsBWs3GPo5tHtKL9VVamW+jrnIp5p4Siu02JibMOvJdozu6t0gv+iz160j9R1lCwDHMWNwf2e2CBhNUF5eXtXQuL0sy3k1fZ3oyRCERq5CW8HB1INokjTsTN5ZLVjYmdnRz6cfal81T3g8gamx6T3uJBhc8n7lyZHM65OmOzwDgz8BKye9VHcuvYDp649zNFlZZj8kwIWPnw7Gy6FhJmFn/+9/pL77HgCOzz+P+9uzRMBoZkTIEIRGqFxbrgSLRCVYZJfe2OfF3tye/j79CfMNo5tHN0yNRLBo9MqKYNdHsO8rQAYbdxi6ENo+qZfqKrUyK/dc4LOos5RVaLE1N2H20HY806Vhei8Asn/4H6nvXQ8YY8fiPmumCBjNkAgZgtDIrDyxkshTkeSU3tjkzcHcgf4+/VH7qnnc43ERLJqSpH1K70XWeaXccTQMnKe33ouEtHzeWB/L8RTl/5/QR1yJ+GsHPBuo9wIg+/sfSJ0zBwDHv43FfaYIGM2VCBmC0MgYSUbklObgaO5If18lWHRp0UUEi6amrBB++wD2LwVksPVQtmN/pHbL5tdURaWWZbsv8MWOBMoqtdhamPDu0EBGdm7ZoL/gs7/7jtT35wLg9MILuM14SwSMZkyEDEFoZIa2GkqgcyBd3LtgYiS+RZukxBil9yI7USk/+jyoPwLLuq0Jcz/xqflMX3+c2EvKInT92roxb0QHWtg37DooWWvXkvbBhwA4jR+P25vTRcBo5sRPMEFoZNys3HCzuvNeL0IjV1oAv70PB5YrZTsveOpLCBigl+rKK7Us+/08X/yWQHmljJ2FCXOGtWfEo14N/ss965u1pH14PWBMmIDb9DdEwBBEyLgTrVZr6CYIeiY+Y6HeXfgdNr8KOcoGgjz2Aqg/ULZm14PTV/N4Y91xTl1RniYc0M6deSOCcLNr+FVcs/7zX9Kur+jrPGkirv/8pwgYAiBCRjVmZsq6AgUFBdjY6GdrY6FxqFqVtLy83MAtEZq80nyIeg8OrVLK9t4w7Evw76eX6soqtCyJPsfineeo0Mo4WJny/rD2DOvoaZBf7Fn/+Q9p8yIAcJ48GdfX/yEChqAjQsZNTExMcHFx4fLlywDY2NhgJJa9fahotVoKCgq4fPkyOTk5aLVamsuCdIIenN8Fm6dCbrJS7jIBwuaCua1eqjt1JZc31sVy+qrSezGwvTsf/CUIN1vD7EGTuWYN6R9/AoDziy/i+o/XRMAQqhEh4xY+Pj6UlpbqgobwcMrJySEtLU03bGJqKp7cEGqhJA80s+FIpFJ28IFhi6FVqF6qK6vQsnhnAkuiz1OhlXG0MmXu8CCGBnsY7Jd65uo1pC6UC54AACAASURBVH9yPWBMCcd12jQRMITbiJBxC0mSaN26NT/88APZ2dm4ubmJb5yHTHl5uS5c5ObmYm1tTYsWLQzcKqHJOLcDNk+DvEtK+fHJMGAOmOtniPXEpVymrz/OmdR8AIZ08OD94e1xsTHcRneZq74mfcECAFxe+jsur74qfk4KdyRCxh0YGRnRtWtXtm3bxtmzZzE1NRXfQA+h8vJyzMzM6NatG46OjoZujtDYFeeA5m04+o1SdlQpvRd+IXqprrSiki9/S2Dp7xeo1Mo4W5sxd3gQQ4I99FJfTWWuWkX6gk8BcHnpJVxefUX8fBTuSmyQdg+pqalcuHCBjIwMKisr9dtAoUEZGRlhZ2eHj48Pfn5+Yu6NcG9nNbBlGuRfUcrdpkD/d8HMWi/VHU/JYfr645xNUyYoP9XRkzlPBeJswN4LgIwVK7j22ecAuLz8Mq6vvmLQ9ggNR2yQpgctWrQQ3eiC0JwVZ8O2WXD8W6Xs1AqGfwW+PfRSXUl5JQt3JLB893m0MrjYmPHhX4IYFGTY3guAjOUruPb59YDxyiu4vvKygVskNAUiZAiCINxJ/K+w5TUoSAUkeOIl6DcbzKz0Ut2R5GymrzvO+WuFAAzv5Mmcp9rjaG2ml/pqI2PpMq4tXAiAy9RXcX3pJQO3SGgqRMgQ/r+9+46Tqrr/P/76bAWW3pa6iKJoVLAbOwY1RI2SqD/zTYwgohJBRVSKhthFRUlUwALSYoqJSrBhwRYLgiWAFDWJwrL0pWzvc35/3FlZCQvLMjP33t338/HYB96ZYc5nj3uH95577jkiUlPxNnhtLCx71jtu1wsunApZJ8aludKKKia9+TXT3/+GiIMOLdK5d+ARnHN4MEZRcx9/nC2PPApAh5E30H7YMJ8rkjBRyBARqbbqZXj5RijaDJYEJ42AM2+F1PjsYPrZmm3c8vdlfJPrjV78/Jiu/O78H9C6mf+jFwBbpk4l99HHAOgwciTth13jc0USNgoZIiJFW2H+aFj+nHfc/hBv9KL78XFprqS8iofe+IoZH36Lc5DZMp37fnYk/Q/LjEt79bFl8hRyJ08GoMOoUbS/+iqfK5IwCmTIMLO5QD/gLefcxXV9TkRkn62cB6/cBEVbvNGLU26AM8ZCanxW0Vz87TZGP7eU1VuLAbjk2G789vwf0KppcBaE2/LYZHKnTAGg48030W7oUJ8rkrAKZMgAHgFmAIP28TkRkbop3AKv3gwr/+EddzgMBk6BrsfGpbni8koefO0rZi9cjXPQqWUTJlx0JGf2Ds6Ou845ch+bTO7UqQB0vOVm2l15pc9VSZgFMmQ45941s377+pyIyF45ByvmegGjeCtYMpx6I5wxGlLisw7Fx99sZfRzy8je5o1e/OL47tx63mG0bBKc0QvnHFsefZStjz8BQMfRo2k35Aqfq5Kw2+eQYWanA7cAxwKdgZ855/6xy2uGR1/TCVgKXOecW7z/5YqI7IfCzfDKKFj1knfc8XBv9KLL0XFprqiskgde+5I5C73t37u0asKEi/pwxiEd4tJefTnn2PLII2x94kkAOo4dQ7vBg/0tShqE+oxkZOAFhxnAC7s+aWaXApOAYcAiYCTwupn1ds5tjr5mSS1tn+OcW1+PmkREauccLH/eG70o2Q5JKXDqKDj9FkiJz50cH/0nl9HPLyNnewkAvzwxi3E/OZQWARq9gGjA+P0f2PrUUwBkjhtL20G6Gi2xsc8hwzk3H5gP1LZe/ShgmnNuZvQ1w4DzgCHA/dH3OKqe9daZmaUDNcc+47P3sogEW8FGeHkUfPWKd5x5JAycCp37xKW5wrJKJry6ij8t8rZ/79q6KQ9c1IdTD24fl/b2h3OOLZMmsXXadAAybx1H28sv97kqaUhiOifDzNLwLqNMqH7MORcxswXASbFsqw7GAbcnuE0RCQrnvAW15o+B0h2QlOqNXJw2CpLjM5rwwb9zGfP8Mtbt8EYvfv3DHoz5yaE0Tw/e9DfnHFsefpit058GIPO222j768t8rkoamlj/5LcHkoFNuzy+CTi0rm8SDSV9gQwzywEucc4t3Ntzu5iAd9mmWgsgp641iEiI5W+Al0fC1695x537eutedDoiLs0VlFZw36ur+MvitQB0b+uNXpx8UPBGL8ALGJsnPsS2GTMAyBz/W9r+6lc+VyUNUfDiNeCcO6s+z+3yujKgrPpYWxGLNALOwZI/w+vjoDTPG73oN9Zb+yJOoxfvfb2Fcc8vY31eKQCDTz6AW37cm4wAjl5ANGA8OJFtM2cCkPm78bT95S99rkoaqlifBblAFbDrsnWZwMYYtyUislPeOm879v+86R13Odobvcj8QXyaK6ng3ldW8rdPvQHSHu2a8cBFffjhge3i0l4sOOfYfP8DbJs9G4BOt/+ONv/3fz5XJQ1ZTEOGc67czD4D+gP/ADCzpOjx5Fi2JSICeKMX//ojvH4blOVDcpq338hJ10FyfEYT3vlyM+Ne+IKN+aWYwRUn9+TmHx9Cs7Rgjl6AFzA2TZjA9jl/BKDTHXfQ5heX+lyVNHT1WSejOdCrxkM9zewoYJtzLhtvHsRsM/sUWIx3C2sGMDMG9YqI7LRjLbx0Pfz3be+463Fw4RToWOcpYPskr7iCu15eyfOfe6MXPdtnMPHiPhx3QNu4tBcrzjk23TeB7X+MBow776TNpf/P56qkMahP7D4OeKfGcfXkytnAYOfcs2bWAbgLbzGuJcAA59yuk0FFROrHOfhsFrwxHsoLIDkdfvRbOGk4JCXHpckFKzdx69wv2FxQhhkMPbUno87uTdO0+LQXK845Nt17H9ufeQaATnffRZtLLvG5KmkszDnndw0JYWYtgby8vDxatmzpdzkiUl/b18CL18G373nH3U/0Ri/aHxyX5nYUl3PnSyuZ+691ABzYwRu9OLZHsEcvIBow7rmX7X/6E5jR+Z67aX3RRX6XJSGUn59Pq1atAFo55/Lr+veCewFRRKSmSAQ+mwFv/A4qiiClKfQfDycOi9voxesrNnLb3OXkFpaRZHDV6Qdy41mH0CQ12KMXAC4SYePdd7PjL3+NBox7aH3Rz/0uSxoZhQwRCb5t33qjF6vf946zTvJGL9odFJ/misq548UVvLjU2+WgV8fmTLy4D0dntYlLe7HmIhE23nUXO/76rBcw7r2X1j//md9lSSOkkCEiwRWJwCfTYMEdUFEMqc2g/+1wwtWQlBSXJud/sYHx85aTW1hOksGwMw7i+v4Hh2L0AqIB48672PFsNGDcdx+tfzbQ77KkkVLIEJFg2vpfb/RizYfecY9T4cLHoO2B8WmusIzfvbiCV5ZtAOCQzOZMvLgvfbu3jkt78eAiETbefgc7/v53MKPL/RNodeGFfpcljZhChogES6QKFj0Jb90FlSWQmgFn3wnHXRm30YtXlnmjF9uKyklOMq7tdxAjftSL9JRwjF5AdcC4nR1/fw6SkryAccEFfpcljZxChogER+5/YN5wWPuxd9zzdLjgMWhzQFya21JQxu/mLWf+cm9B4kM7teChS/pyRNdWcWkvXlwkwobx48l7/gUvYDzwAK1+er7fZYkoZIhIAESq4OOp8PY9UFkKac3hnLvh2CsgDvsOOed4cel67nhxBduLK0hJMoaf2YvhZ/YiLSU+oyXx4iIRNvx2PHkvRAPGgw/S6vzz/C5LBFDIEBG/bfka5l0LOZ94xweeCRc8Cq2z4tLc5oJSfjt3OW+s9NYHPKxzSx66pA+HdwnX6AWAq6ryAsbcuV7AmPggrc5TwJDgUMgQEX9UVcLCyfDOfVBVBukt4Zx74JjL4zZ6MW/Jem5/cQV5JRWkJhsjzjyYa888iNTkcI1eQDRg3HobefPmQXIyXSc+SMtzz/W7LJHvUcgQkcTbvAr+cS2s/9w77nUW/PQRaNUtLs1tyi/ltrlfsGDVZgCO6NqSiRf35bDO4Vz91wsYt5I370UvYDz8EC0HDPC7LJH/oZAhIolTVQkfPQLv3g9V5ZDeCgbcB0f9Km6jF89/vo67XlpBfmklqcnGyLMO4erTDwzl6AV4AWP92HHkv/RSNGA8TMsBP/a7LJHdUsgQkcTYtMIbvdiwxDs++Mfw0z9Ayy5xaW5jXinjXljGO19tAaBPt1ZMvLgvvTu1iEt7ieAqK72A8fLLkJLiBYwfn+N3WSK1UsgQkfiqqoAPfg/vPQiRCmjSCn7yIPS5NG6jF3//NIe7X1lJQWklaclJ3Hj2IVx1Wk9SQjp6AdGAMWYs+a+84gWM30+i5dln+12WyB4pZIhI/Gz8whu92LjMO+59Lpz/e2jRKS7Nrd9RwtgXvuCfX3ujF0d1b83Ei/twcGZ4Ry8gGjBGjyH/1VchJYVuf/g9Lc46y++yRPZKIUNEYq+yHN5/GN5/CCKV0LQN/GQiHHlx3EYv/vrJWu59ZRWFZZWkpSRx8zmHcOWpB5KcFPv2EslVVrLullsomP8apKZ6AaN/f7/LEqkThQwRia0NS73Ri03LvePDfgrnPgwtMuPSXM72Ysa98AXv/zsXgGOyWvPgxX3p1bF5XNpLJFdRwbpbRlPwWjRgPPIILX50pt9lidSZQoaIxEZlGfxzIrw/CVwVNGsH5z4Eh/8sLqMXkYjjz4uzmfDqKorKq0hPSeKWH/fmilN6hn70AqIB46abKXjjDSw1la6PPkKLMxUwJFwUMkRk/6373NtzZPNK7/gHA72A0bxDXJpbu62YMc8v46P/bgXg+APa8ODFfenZPiMu7SWaq6hg3aibKHjzTS9gPPYoLfr187sskX2mkCEi9VdZ5q158eEj0dGL9nDew3D4wLg0F4k4nlm0hvvnf0lxeRVNUpMYM+BQBp10AEkNYPQCwJWXs+6mmyh4cwGWmkq3yY/R/Iwz/C5LpF4UMkSkfnI+8/Yc2fKld3zERd7kzox2cWluzdYiRj+3jEXfbgPghJ5tefCiPhzQQEYvwAsYOaNGUbjgLSwtjW5TJtP8tNP8Lkuk3hQyRGTfVJTCu/fBR4+Bi0BGRzh/kjfBM05KK6q4cMqH7CiuoFlaMmN/ciiXndijwYxeQDRgjLyRwrffjgaMKTQ/7VS/yxLZLwoZIlJ3axd7cy9yv/aO+1wKA+6HZm3j2mxeSQU7iitIMnjthtPJatcsru0lWqS8nHU3jKTwnXew9HQvYJx6it9liew3hQwR2buKEnj7Hlg4BXDQvJO3qNahidn1s7wyAkBaSlLDDBjX30Dhu+96AWPqFJqfooAhDYNChojs2ZqF3ujFtv96x31/6W1q1rRNwkqojDgAUpPCuyz47kTKy1l33fUUvvcelp5O98enknHyyX6XJRIzChkisnvlxfD23fDx44CDFl28Dc0OSfyOn5VV3khGSnLDmYMRKSsj5/rrKXrvn1iTJl7AOOkkv8sSiSmFDBH5X6s/9EYvtn/rHR99GZxzLzRt7Us5FVXRkYwQb3BWU6SsjJwR11H0/vtewHjicTJ++EO/yxKJOYUMEdmprBDeuhMWP+Udt+wKP30UDvZ3M66K6EhGQwgZkbIycoaPoOiDD7CmTen+xBNknHiC32WJxIVChoh4vv0nzBsBO9Z4x8cOhrPvhiYtfS0LoDLSMC6XREpLvYDx4YdewHjyCTJOUMCQhkshQ6SxKyuAN2+HT5/2jlt1hwsehYN+5G9dNVRfLkkJ8boYkdJScq4dTtFHH2HNmpH15BM0O/54v8sSiatAhgwzmwv0A95yzl28y3OrgXwgAmx3zmnHIJH6+uZdmHcd5GV7x8cNgbPvgvQWvpa1q8qQz8mIlJSQM3w4RR8t9ALGU0/S7Ljj/C5LJO4CGTKAR4AZwKBanj/ZOVeYwHpEGpbSfHhzPHw2yztunQUXTIYDg7lHRpjnZERKSlj7m2sp/vhjkpo1o/u0p2h27LF+lyWSEIEMGc65d82sn991iDRI/3kLXrwe8nO84xOuhv63Q3pzf+vag4qQ3sIaKS72AsaiRV7AmD6NZscc43dZIgmzz78WmNnpZvaSma03M2dm/7PdopkNN7PVZlZqZovMLJYzmxzwnpl9Yma/iuH7ijRspXnexM5nfu4FjDYHwOBX4NyJgQ4YEM7FuCLFxawd9hsvYGRk0H36dAUMaXTqM5KRASzFu5zxwq5PmtmlwCRgGLAIGAm8bma9nXObo69ZUkvb5zjn1u+l/VOdc+vMrDOwwMy+cM4tq8f3IdJ4fP0GvHQDFKwHDE4cBv3HQ1o4djD97nJJSjhGMiLFxay9ZhjFn3wSDRjTaHb00X6XJZJw+xwynHPzgfkAZrs94UcB05xzM6OvGQacBwwB7o++x1H1rBfn3LronxvM7FXgGOB/QoaZpQPpNR4K1kw2kUQo2Q6v3QpL/+wdtz0ILpwCPcK1suTOu0uCP5IRKSryAsann5LUvDlZ06fR9Kh6f+SJhFpMz1gzSwOOBRZUP+aci0SP9/tTzcwyzKxF9L+bAz8CVtTy8nFAXo2vnP1tXyRUvpoPU34YDRgGJ42AYR+ELmDAzmXFUwM+J6OqsIjsq6/ZGTCenq6AIY1arCd+tgeSgU27PL4JOLSub2JmC4C+QIaZ5QCXOOcWApnA3OgISjLeiMkntbzNBLzLNtVaoKAhjUHxNnhtLCx71jtud7A3epF1or917YeKSPBHMqoKi1h79dWUfP45SS1aeAGjTx+/yxLxVVDvLtntGsbOuW/wwkdd3qMMKKs+ruXSjkjDsuplePlGKNoMluSNXpx5K6Q29buy/fLdSEZKMENGVWEha6+6mpJ//csLGDOepumRR/pdlojvYh0ycoEqvBGHmjKBjTFuS0SqFW2F+aNh+XPecfveMHAqdGsYCz59N/EzgCt+VhUWsnboVZQsWUJSy5ZkPf00TY88wu+yRAIhpr8WOOfKgc+A/tWPmVlS9HhhLNsSkaiV82DqiV7AsCQ4dRRc888GEzCgxsTPgM3JqCooYO2VQ72A0aoVWTNmKGCI1LDPIxnRCZe9ajzU08yOArY557Lx5kHMNrNPgcV4t7BmADNjUK+IVCvKhVdvhhVzveMOh8HAKdC14a0mWfldyAjO5ZKqggKyhw6ldOmyaMB4mqaHH+53WSKBUp/LJccB79Q4rp5cORsY7Jx71sw6AHcBnYAlwADn3K6TQUWkPpzzgsWrN0PxVrBkOG0UnH4LpKTv/e+HUPUurGkBCRlV+flkD72K0mXLSG7ViqyZM2jygx/4XZZI4NRnnYx3gT2OWTrnJgOT61mTiNSmcDO8chOsetE7zjzCu3OkS8O+TbK8elnxAMzJqMrPJ/vKoZR+8YUXMGbNpMlhh/ldlkggBfLuEhHZhXOw/Hl49RYo2QZJKXDazXDaTZCS5nd1cReUyyVVeXlewFi+nOTWrb2AcWid784XaXQUMkSCrmAjvDwKvnrFO+50JFw4FTo3njUYgrAYV1VeHtlDrqR0xQqS27TxAkbv3r7VIxIGChkiQeWct6DW/DFQugOSUuGM0XDqjZCc6nd1CVW9GJdfW71X7djhBYyVK6MBYxZNeh/iSy0iYaKQIRJE+Rvg5ZHw9Wvecee+MPBxyGycdy9UVPq31XvVjh2sGTKEspWrSG7b1hvBOEQBQ6QuFDJEgmbpX72FtUrzIDkNzhgDp9zQ6EYvavJrq/fK7dvJHnIlZatWkdyuHT1mzST94IMTWoNImClkiATN5lVewOhyjLdqZ0fduVC94mciRzIqt28n+4ohlH35Jcnt23sBo1evvf9FEfmOQoZI0PQbB627wzGDIVmnKOy8uyRRczIqt28ne/AVlH31lRcwZs8i/aCDEtK2SEOiTzCRoEltAscP9buKQKlI4N0lldu2eQHj669J7tCeHrNnk37ggXFvV6QhCsbyeSIie5Cord4rt24le9BgBQyRGNFIhogEXmUC5mRUbt1K9uDBlP37P6R06EDW7NmkH9gzbu2JNAYayRCRwKuekxGvvUsqc3NZM2iQFzA6diRrjgKGSCwoZIhI4H23d0kcQkblli2sGTSY8v/8l5TMTHrMmU16TwUMkVjQ5RIRCbzqXVhjfbmkYvNmsgdfQfk335DSqRM9Zs8irUePmLYh0pgpZIhI4H13C2sMJ35WbN5M9qDBlH/7rRcw5swmLSsrZu8vIrpcIiIhEOtbWCs2bSb78kFewOjcWQFDJE4UMkQk8CpiuNV7xaZNZA8aRPnq1aR0UcAQiSddLhGRwIvVVu8VmzZ5Ixhr1kQDxhzSunWLRYkishsKGSISeLFYjKti40bWDBpExZpsUrt0IWvOHNK6dY1ViSKyG7pcIiKBVz2SkZZSv5GMig0bWHN5NGB07aqAIZIgChkiEnjfzcmox0hGxfr1XsDIzia1WzdvDoYChkhCKGSISODVd6v3inXrvICxdi2p3bvTY85sUrsqYIgkiuZkiEjgVUb2fav37wLGunU7A0bnzvEqUUR2QyMZIhJokYijah9DRnlOjYCRlUWPP85RwBDxgUYyRCTQKqJLikPdLpeU56wj+/LLqVi/ntQeWfSYM4fUzMx4ligitVDIEJFAq15SHPa+rHh5Tg5rLr+cyvUbSOvRg6w5sxUwRHykyyUiEmg1Q8aeRjLK165lza+jAeOAA8jSCIaI7xQyRCTQvne5JGn3IaM8O9sLGBs2kNazZ3QEo2OiShSRWihkiEig1dwczex/Q0b5mjWsuXwQlRs3knbggWTNnkVqRwUMkSDQnAwRCbTKPSzEVb56NWsGDaZy0ybSDjqIHrNmktKhQ6JLFJFaKGSISKDVthBX+erV3gjG5s2k9TqIHrNmkdK+vR8likgtAnm5xMzmmtl2M3tuN8/dbGYrzGy5mV3mR30ikjjVC3Gl1Vgjo+ybb705GJs3k35wLwUMkYAKZMgAHgEu3/VBMzsS+CVwLHA8MMLMWie4NhFJoPLK749klH3zLdmDBlG5ZQvpBx9MlgKGSGAFMmQ4594FCnbz1GHAQudcqXOuBFgKDEhkbSKSWJU1tnkv++Yb1gy63AsYhxxC1uxZpLRr53OFIlKbfQ4ZZna6mb1kZuvNzJnZwN28ZriZrTazUjNbZGYnxKZclgP9zKy1mbUB+gHa7UikAave5r17wUbWXD6Iqi25pPfuTdasmaS0betzdSKyJ/WZ+JmBN4IwA3hh1yfN7FJgEjAMWASMBF43s97Ouc3R1yyppe1znHPra2vYObfSzB4F3gbygI+Bqnp8DyISEhVVjqz8jdy48CmqSvJJP/RQsmbOIKVNG79LE5G92OeQ4ZybD8wHdnvPOjAKmOacmxl9zTDgPGAIcH/0PY6qZ704554Enoy+93Tg37t7nZmlA+k1HmpR3zZFxD+Rb//L/R8+QcuyQtIPO4ysGU8rYIiEREznZJhZGt6kzAXVjznnItHjk2LURsfon72BE4DXa3npOLzRjuqvnFi0LyKJU/r117S67QbalBWyrkOWAoZIyMR64md7IBnYtMvjm4BOdX0TM1sA/B0418xyzKxmQJlnZiuBZ4ArnHOVtbzNBKBVja9udW1fRPxX+tXXZA++guS8Hfy7VVdm/myUAoZIyARyMS7n3Fl7eK5OIyLOuTKgrPq4lks7IhJApV99RfbgK6javp2yAw/h1kMv4/CMln6XJSL7KNYjGbl4EzF33fowE9gY47ZEpAEq/fJLsgcNpmr7dpoccQTZ4x6gMK0ZqcmBvONeRPYgpmetc64c+AzoX/2YmSVFjxfGsi0RaXhKV63yRjB27KDJkUeSNeNpyppmAHve5l1EgmmfL5eYWXOgV42HeprZUcA251w23u2rs83sU2Ax3i2sGcDMGNQrIg1U6cqVZF8xhKq8PJr06UPW9Gkkt2xJZWQHsPsN0kQk2OozJ+M44J0ax5Oif84GBjvnnjWzDsBdeJM9lwADnHO7TgYVEQGgZMUKsodcSSQvjyZ9+5A1fTrJLby7zqsX40pL0UiGSNjUZ52Md4E9nu3OucnA5HrWJCKNSMnyFWRf6QWMpn370n36tO8CBkD5HrZ6F5FgC+TdJSLSOJQsX0H2kCFE8vNpetRRXsBo3vx7r6msZat3EQk+/WogIr4o+eKLnQHj6KPpPn36/wQM2LlBWqpGMkRCR2etiCRcybJl3hyM/HyaHnMM3adNI7l5xm5fWxEdyUjVnAyR0FHIEJGEKlm61AsYBQU0Pe5Yuj/1VK0BA3aGDM3JEAkfnbUikjAlS5aQfeVQIoWFNDvuOLKefHKPAQOgMjrxM1VzMkRCRyFDRBKi+F//2hkwjj+e7k8+QVLGngMGeFu9A6RoxU+R0NFZKyJxV/z5v1g79CoiRUU0O+GEOgcMgMpIdE6GQoZI6OisFZG4Kv78c9YOHeoFjBNPpPsTj5PUrFmd//53Ez+TdLlEJGwUMkQkboo/+8wbwSguptkPf7jPAQN0uUQkzHTWikhcFH/6KdlXXU2kuJiMk0+i++NTSWradJ/fp3oxLk38FAkfhQwRibniTz4h++prcMXFZJx8Mt2m1i9gAFRUL8alkQyR0NFZKyIxVbR48c6AccopdJs6haQmTer9fhWVWlZcJKwUMkQkZooWLWbtNcNwJSVknHoq3aZM3q+AAVpWXCTMdNaKSEwUffwxa6+5xgsYp50Wk4ABNVb81EiGSOgoZIjIfitauJC1w36DKy0l44zT6Tb5MZLS02Py3jtX/NTHlUjYaKt3EdkvRR99xNrfXIsrK6P5GWfQ9bFHSUpLi9n7V+juEpHQ0q8GIlJvhR9+uDNg9OsX84ABO+8u0QZpIuGjs1ZE6qXwgw/JqQ4YZ55J10cfiXnAgJ3rZGhOhkj4KGSIyD4rfP8Dcq69FldeTvP+/en2yB/iEjBg55yMNM3JEAkdnbUisk8K//lPcoYP9wLGWf3p9vtJWJwCBtS8u0QfVyJho7NWROqs8L33yBk+AldeTouzz6LbpPgGDICKiC6XiISVQoaIZmVdAwAADLtJREFU1EnBu++SM+I6XEUFLc4+m64JCBhQ4xZWTfwUCR2dtSKyVwXvvEPOddd7AeOcc+g66WEsNTUhbVfvwpqaopEMkbBRyBCRPSp4+x1yrr8BKipoMWAAXR9+KGEBA2rMydBIhkjo6KwVkVoVvPUWOTdEA8ZPBtD1oYkJDRigrd5FwkwhQ0R2q2DBAnJuGAkVFbQ891y6TpyIpSR+keDvFuPS3SUioaOzVkT+R/6bb5Iz8kaorKTleefR5cEHfAkYoJEMkTBTyBCR78l//Q3W3TjKCxjnn0+XB+73LWBURRzRgQzdXSISQjprReQ7+a+9zrpR0YBxwU99DRiwc9InaJ0MkTAKXMgws+5m9q6ZrTSzZWZ2SY3nWpvZp2a2xMyWm9lVftYq0pDkv/Ya6266CaqqaHXhBXSZMAFLTva1psrqYQy01btIGAVxq/dKYKRzbomZdQI+M7NXnXNFQAFwunOu2MwygOVm9oJzbquvFYuEXP6rr7LultHRgHEhne+71/eAATvnY4BChkgYBS5kOOc2ABui/73RzHKBtkCRc64KKI6+NB2w6JeI1FPeK6+wfvQYL2AMHEjne+8JRMAAKI+GDDNITtKpLhI2+/yrgZmdbmYvmdl6M3NmNnA3rxluZqvNrNTMFpnZCfUpzsyOBZKdc2trPNbazJYCOcBE51xufd5bRCDvpZdZXz2C8fOfBypggJYUFwm7+py5GcBSYPjunjSzS4FJwJ3AMdHXvm5mHWu8pnpOxa5fXWq8pi0wB7i65vs753Y45/oCPYFfmllmPb4HkUYv76WXWD9mDEQitLr4Ijrfc3egAgbsDBma9CkSTvt8ucQ5Nx+YD2C22xN/FDDNOTcz+pphwHnAEOD+6Hsctac2zCwd+Adwv3Puo1rq2BQd0TgNeK6W90iv8VCLPX5jIo1I3osvsn7sOIhEaH3JxXS6804sgKMF1Tuwaj6GSDjF9Mw1szTgWGBB9WPOuUj0+KQ6vocBs4C3nXN/3OW5TDNrEf3vVsDpwFe1vNU4IK/GV86+fC8iDVXevHmsHzM2GjAuCWzAgJ23sGohLpFwivUnS3sgGdi0y+ObgE51fI9TgEuBgdHLKkvM7Mjocz2A96MjGO8DjznnvqjlfSYArWp8dav7tyHSMO2Y+w9vBMM5Wl96KZ3uvCOwAQNqXC4JcI0iUrsg3l3yAbWEH+fcYmCPl1pqvLYMKKs+ruXSjkijseOFuWy47TYvYPzfL+g0fnygAwbU2IFVIxkioRTrkJELVAG7TsbMBDbGuC0RqaMdzz/Pht+OB+do88v/I3P8+FAE7+rFuNI0J0MklGJ65jrnyoHPgP7Vj5lZUvR4YSzbEpG62fHcczsDxq9+FZqAAVBRqZEMkTDb55EMM2sO9KrxUE8zOwrY5pzLxrt9dbaZfQosBkbi3fY6Mwb1isg+2P63v7Hxd7cD0Oayy8i87dbQBAyosc17wC/riMju1edyyXHAOzWOJ0X/nA0Mds49a2YdgLvwJnsuAQY453adDCoicbT92b+x8fZowLj812SOGxeqgAHa5l0k7OqzTsa77GUpb+fcZGByPWsSkf20/a/PsvGOOwBoO+hyOo4dG7qAAVBRveKn5mSIhFLg7i4Rkf2z/S9/YeOddwHQdtAgOo4dE8qAAbq7RCTsFDJEGpBtf/4zm+66G4C2V1xBx9G3hDZgAFRqxU+RUFPIEGkgtj3zJzbdcw8Aba8cQsebbw51wICdl0tStAOrSCgpZIg0ANv++Ayb7r0XgHZDr6TDTTeFPmBAjV1YNZIhEkoKGSIht23OHDbdNwGAdlddRYdRNzaIgAE19y5RyBAJI4UMkRDbOmsWm+9/AIB211xDh5E3NJiAAZr4KRJ2ChkiIbV15iw2PxANGMOuocMNDStgwM5lxbUYl0g4KWSIhNDWGTPZ/OCDALS/9je0v+66BhcwYOdiXGkpDe97E2kMFDJEQmbr00+zeeJDALS/9lraXzeiQQYMgHJt9S4SagoZIiGSO20aWx72VvJvP2IEHUYM97mi+KrUnAyRUFPIEAmJ3KemsWVSNGBcN4IOwxt2wICdczJ0d4lIOClkiIRA7hNPsuUPfwCg/fXX0eHaa32uKDEqtEGaSKgpZIgEXO4TT7DlD48A0GHkDbQfNsznihLnu1tYNSdDJJQUMkQCbMvUqeQ++hgAHW68kfbXXO1zRYm1c8VPjWSIhJFChkhAbZk8hdzJkwHoMGoU7a++yueKEu+7vUs0J0MklBQyRAJoy2OTyZ0yBYCON99Eu6FDfa7IH9qFVSTcFDJEAuZ7AeOWW2h35RCfK/KPJn6KhJtChkjApHRoD0DH0aNpN+QKn6vxV4UW4xIJNYUMkYBp84tf0PToo2nSu7ffpfhOi3GJhJt+PRAJIAUMT/ViXGmakyESSjpzRSSwyis1kiESZgoZIhJY2updJNx05opIYFXq7hKRUFPIEJHAqqjSBmkiYaYzV0QCq0J3l4iEmkKGiASWtnoXCTeduSISWDt3YdVIhkgYKWSISGBVak6GSKjpzBWRwNq5d4k+qkTCSGeuiASWJn6KhFuj27skPz/f7xJEpI7KSgqJlEcoLSokPz3idzkijVZ9/+0051yMSwkmM+sK5Phdh4iISIh1c86tq+uLG1PIMKALUOB3LXXUAi8UdSM8Ncea+sCjflAfVFM/qA+q+dEPLYD1bh+CQ6O5XBLtlDqnL795mQiAAudco7zGoz7wqB/UB9XUD+qDaj71wz63o4mfIiIiEhcKGSIiIhIXChnBVQbcGf2zsVIfeNQP6oNq6gf1QbVQ9EOjmfgpIiIiiaWRDBEREYkLhQwRERGJC4UMERERiQuFDBEREYkLhYyQMrO5ZrbdzJ7b5fHWZvapmS0xs+VmdpVfNSZCbf2wt+caKjO72cxWRP/fX+Z3PX4ws97Rn//qrxIzG+h3XYlmZqvNbFm0D97xux4/NLbPw9r4+Vmou0tCysz64S3xOsg5d3GNx5OBdOdcsZllAMuB45xzW/2pNL5q64e9PdcQmdmRwGzgZMCAd4ABzrkdvhbmIzNrDqwGejjninwuJ6HMbDVwhHOu0O9a/NLYPg9r4+dnoUYyQso59y67Wa/eOVflnCuOHqbj/WPTYPfJrq0f9vZcA3UYsNA5V+qcKwGWAgN8rslvFwBvNbaAIZ7G9nlYGz8/CxUy4sDMTjezl8xsvZm53Q3Vmtnw6HBmqZktMrMTYth+azNbird5zkTnXG6s3nsf6/C1H4ImAf2xHOgX/f/fBugHdI1R+TGT4J+L/wc8u38Vx16C+sAB75nZJ2b2q9hUHluJ6IegfB7WpqF/TjaaDdISLAPvt8gZwAu7PmlmlwKTgGHAImAk8LqZ9XbObY6+Zgm7//9zjnNu/Z4ajw6P9zWzTOAFM3vOObdpf76hevK1HwIo3v2x0sweBd4G8oCPgap4fCP7KSE/F2bWEu/S0S/i8U3sp0T0wanOuXVm1hlYYGZfOOeWxefbqbe490OAPg9r07A/J51z+orjF95vEwN3eWwRMLnGcRLeDrFj9/G9+wHP7eU1U4GLG2s/1KWPGlp/1Pj704Hz/P5effy5+DXwjN/fY0B+FiYCg/3+XgPQD4H4PPSjD/z6LNTlkgQzszTgWGBB9WPOuUj0+KQYvH+mmbWI/ncr4HTgq/1931iLdz+ETaz6w8w6Rv/sDZwAvB7bSuMrxj8XgbxUsjex6AMzy6jxOdAc+BGwIvbVxk+M+iEUn4e1aQifk7pcknjtgWRg1+G6TcChdX0TM1sA9AUyzCwHuMQ5txDoATxlZtUTnB5zzn0Rk8pjK979sMfnAigm/QHMi36YFgFXOOcqY1RfosTq56IVXsi6KHalJUws+iATmOt9DJAMTHPOfRKzChMjFv0Qls/D2sT9czLeFDJCyjl3Vi2PLwaOSnA5vqmtH/b2XEPlnAvFbzfx5pzLw/uHtlFyzn2D949Ko9bYPg9r4+dnoS6XJF4u3mS8XT8AM4GNiS/HN+qH71N/eNQP6oNq6ocG0AcKGQnmnCsHPgP6Vz9mZknR46AO5cec+uH71B8e9YP6oJr6oWH0gS6XxEF0olWvGg/1NLOjgG3OuWy825Fmm9mnwGK8W5IygJkJLzaO1A/fp/7wqB/UB9XUD42gD/y+ZachfuHdKuR28zWrxmtGAGuAMrxblE70u271g/pD/aA+UD+oD2L5pb1LREREJC40J0NERETiQiFDRERE4kIhQ0REROJCIUNERETiQiFDRERE4kIhQ0REROJCIUNERETiQiFDRERE4kIhQ0REROJCIUNERETiQiFDRERE4kIhQ0REROJCIUNERETi4v8Dp69z/T7aPz4AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "print(\"Solution Found:\",x0) # Print the solution\n", "print(\"History:\") # Print the values of the function throughout the optimization\n", "for i in range(0,vals.size):\n", " print(i,\" \",vals[i])\n", "\n", "dis = np.linspace(1,Niter,Niter+1) # Error analysis \n", "errors = abs(vals-analytic) # Compute differences between current points and the optimum\n", "if(dis.size>10):\n", " dis2 = dis[0:10]\n", "else:\n", " dis2 = dis\n", "sq = 0.1**dis2 # Construct curve of order 1\n", "sq2 = 100*sq**2 # Construct curve of order 2\n", "sq3 = 100*sq**3 # Construct curve of order 3\n", "\n", "\n", "plt.figure(2)\n", "plt.loglog(errors[:-1:],errors[1:],label='Errors') # Plot the errors in log-log plot\n", "plt.loglog(sq,sq2,label='Order 2') # Plot order curves for comparison\n", "plt.loglog(sq,sq,label='Order 1')\n", "plt.loglog(sq,sq3,label='Order 3')\n", "plt.legend(loc='best', shadow=True, fontsize='large') # Show legend\n", "plt.show() # Show plot" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "celltoolbar": "Raw Cell Format", "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.6.7" } }, "nbformat": 4, "nbformat_minor": 2 }