{ "cells": [ { "cell_type": "code", "execution_count": 93, "metadata": {}, "outputs": [], "source": [ "# Import libraries\n", "import time\n", "import numpy as np\n", "\n", "import numpy.polynomial.polynomial as poly\n", "from scipy.interpolate import lagrange\n", "\n", "import matplotlib.pyplot as plt\n", "plt.rcParams['text.usetex'] = True" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## MA934\n", "\n", "## Interpolation\n", "\n", "### Lagrange interpolation (and conditioning)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "The following code provides the building blocks associated with the third set of notes (on Lagrange interpolation). It features two approaches:\n", "* a direct implementation of the strategy in (3.1) which introduces the relevant polynomials\n", "* usage of the scipy functionality to achieve the same goal in an in-built manner\n", "\n", "The functions we test below are the Gaussian and Runge functions found in (3.2) and (3.5), respectively." ] }, { "cell_type": "code", "execution_count": 96, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvoAAAH6CAYAAAB29+hIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABzg0lEQVR4nO3deXxU1f3/8ddJCJCQwITFoHVjYkVwTwJqV5dJS4uIrQmIVqq2JLa12pWU77fa1i6YtH7bWm2b0GKl/amQWBcU0UTFbgpZ6opazeBaQTQMJCRsyfn9MTMxyySZJJPcyc37+XjMA3LXzz33zp3PnDnnXGOtRURERERE3CXB6QBERERERCT2lOiLiIiIiLiQEn0RERERERdSoi8iIiIi4kJK9EVEREREXEiJvoiIiIiIC41xOoBYMcZ4gcIOk7xAkbXWH+X6HmAFUB+alGmtLYppkCIiIiIiw8QViX4oyc/rmJgbY/KAWmNMdpTJfjlQGF7WGOM1xlRaa3OHJmoRERERkaHjlqY7hV0nWGsrAE+keV2FvhT4O34h6JDw58UuTBERERGR4eGWRB9gcYRpAYLJfjTr1kaYXkkUXxREREREROKNKxJ9a22RtTa747RQm3sPwSY5ffEBkZr3+IGcwcYnIiIiIjLcXJHo92AVUGatreptoQ5fCBoizA4Q3S8CIiIiIiJxxRWdccOMMVkEa+fnAGtD7fT7MjmK7XqstYFBhiciIiIiMmxclehba+uAutAoPEXGmMnW2rKh2JcxpgAoABg/fnz20UcfPRS7GXXa2tpISHDzD03DS+UZWyrP2FJ5xo7KMrZUnrGl8oyd//znP+9Za6dFu7yrEv2w0Ig5hcaYXaHa+JIh2EcZUAYwc+ZM+/LLL8d6F6PSpk2bOPvss50OwzVUnrGl8owtlWfsqCxjS+UZWyrP2DHGvN6f5d3+9aoMKO5jmUht8ztRsx0RERERGWlGfKJvjPEYYypD7fO7ej+0jLen9UNJfIDgk3S78obmiYiIiIiMKCM+0SeYjPuInKhPCf3bV619DZE75WYCvY7aIyIiIiISj0Z8oh/qgFvSwwg7PqAqiqY35UBuD+uvHVyEIiIiIiLDb8Qn+iHVoVFw2hljwrX8hV2m1xtjOj0FN9Sx1tuxiU+oKVBDlEN0ioiIiIjEFVeMumOtrTDGZBljigm2y59CMMmfEaE2P0DkpjznASuMMfWhvzOttZFq+UVERERE4p4rEn34YAz9KJbL7mF6ACiKcVgiIiIiIo5wS9MdERERERHpQIm+iIiIiIgLuabpzkiyf/9+GhoaaGxspLW11elw4sqkSZN48cUXnQ7DNVSesaXyjK14L8+EhATGjx9Pamoq6enpJCSobkxERhYl+sNs//79vPHGG6Snp3PssceSlJSEMcbpsOJGY2MjaWlpTofhGirP2FJ5xlY8l6e1lra2NpqbmwkEAuzZs4ejjjqKMWP0sSkiI4eqJ4ZZQ0MD6enpTJ06lbFjxyrJFxGJQ8YYEhMTSUtL48gjj2TcuHE0NPT17EURkfiiRH+YNTY2MnHiRKfDEBGRKBljmDJlCrt373Y6FBGRflGiP8xaW1tJSkpyOgwREemHsWPHcujQIafDEBHpFyX6DlBzHRGRkUX3bREZiZToi4iIiIi4kBJ9EREREREXUqIvIiIiIuJCSvRFRERERFxIib7EhUAggDGGhQsX9rhMRUUFxhhKSkqGMTL3yM3NJTMzM6bbDJ+33NzcmG53NCspKcEYQ1VVldOhDIuysjKMMdTV1TkdioiI6yjRFxGJIx6PB6/X63QYIiLiAkr0RUTiSEFBAfX19fh8vgFvo6ysbNT8ItAfKhcRGW2U6IuIuExRURHl5eVOhxF3VC4iMtoo0RcRERERcSEl+uIKgUCAwsJCMjMzSU9PJz8/H7/f32mZiooK0tPTqauro7CwkPT0dNLT0yksLATA7/eTm5uLMYbMzEzKysoi7ic/P5/09HQyMzMpKiqKOsbB7j8QCFBUVERmZmb7Mj3tv6qqiuzsbI4++miys7OpqKigoaGhx7Ib6DH1R3/ir6urIzs7G2MM2dnZlJWVkZ+fT2ZmZnvH33B5+v1+qqqqyM3NbS/LaPfV8ZyElw9fP5E8/fTTfcYV3v9AyzRS59Ro48zPz8cYQyAQaN9O1w7sfcXWW7mGr6u6ujpKSkrayyI3N5dAINDtWKJ5X0YS6/MXTbmIiLiREv04E/4AGmmvWNm2bRtVVVURX9XV1RHX8fv9zJgxg6qqKgoLCykuLsbv95OZmdmpPW5DQwOBQIDzzjsPgOLiYnJyctqTtdzcXHJzcykuLqahoYHCwsJOyVZ4P3V1dRQXF1NYWEhJSUl7EtSXwe5/3bp1VFVVkZeXR3l5OXl5eRH3X1dXR25uLnV1dVx++eUUFhaycuXKiKOaDPaY+iPa+P1+P9nZ2Xi9Xmpra8nJyaGwsBCv10tpaWmnhC8QCFBaWkpubi4NDQ3tyXa0+wqfk/z8/PbkMicnh4qKim6Jpd/v5xOf+ESfcQ1FmUYbZ3FxMZWVlQD4fD4qKyuprKykoKCgX7H1VK6BQKD9S1h1dTWLFy8mKyur/QtA1/KK5n0ZSazPX1/lIiLiWtZavQb5Ov744220tm7d2ut8YES+BmvXrl1R76u4uLjTuj6fz3o8nm7bzMrKsl6vt/3v0tJSC9i8vLxOy3k8HgvY8vLy9mmVlZXd9uXz+Sxgd+3a1W2bHaf1ZLD7jyQcU9fjBmx9fb3ds2dP+3Sv19upPGJxTOHz5vP5+lw22vgLCgq6nU+Px9Ot3MrLy9uviY5l1599hY+1a/xAt7KKNq5YXSe1tbUDijM8vaCgoNv0aGLrrVzD80pLSyNut7KystO0vt6Xe/bsiXi8kQz2/IWnRyqXaPV1/3bS448/7nQIrqLyjC2VZ+wANbYfOapq9ONMf05ePL1i5ZxzzulxH5E60YWbF6xYsaLbvHANYkVFRafpixcv7vR3eCjDvLy89mk5OTkAvP/++0CwJjNcw9jQ0IDf78fv97ev25+RPAay/56Ea1rDzSHCNa4FBQXdhmj0eDyd/o7lMQ1U1/jD/+8ae05OTo9NPvLy8jqVXX/21XVeWFZWVremTn6/n2OPPbbXuIa6TKOJsyf9ja23cu16fkpLSwHa36MDeV/2ZbDnT0RkNBrjdAAigxFujpKVldVtXjhZrq6u7pSwdE14J0+e3G1aVzU1NUCwTXCkBMXv90dsGhNpTPSB7L/jfioqKqisrGxP1LrOB6J6MFY0xxRrfcUPwSQyHFvHWHsabrLrF6f+7Css0vUzkLiGukyjibMn/Y2tp3KNJHyNh7cR7fvy05/+dI/bjPX5ExEZjVSjL64QqSNgTzV6kydPHvB+KisrI/7asHz5cpYtW0Z2dnanV6QOhAPdf0VFBZmZmaxdu5b8/HzKy8tZvnx5xGWj/eLQ1zHFUrTxh9tbh/soFBYWEggEItYOQ/fa5f7sKyyac9KfuIaqTAdz7fY3tlg8tKs/78uOhuL8iYiMRkr0ZUQL1+RF6qgbrlWcM2fOoPcTroWMVGsfVltbG1Vzo4FatmwZeXl51NbWUlBQQFZWFlOmTOm0TDg5q62t7XN70RxTLEUTfzgen89HTU0N2dnZ1NTUUFlZ2a9a22j31R91dXWcc845vcY13GXam64J9VDGFq5tD5fFYN+XQ3H+wtSkR0RGEyX6MqJ5vV6ysrIoKyvrVntYVFSEx+OJqv12XzweDz6fj5UrV3bbTyAQiFhzGUvhfXStZe2aSHk8nojlEQgEujV9GM5jijb+8LTc3Nz2L061tbX9ekpsf/bVH9XV1Zxzzjm9xuX0ddIxjq77imVsXa+l8Gg44eY+g3lfDtX5g8jlIiLiZmqjLyNeeXk52dnZzJgxo70ZRWlpKX6/v31IvVgoLS3ttB+Px0NtbS1lZWXtQwAOlXBb/5KSEgKBALm5uaxduzZiW+vi4mJyc3OZMWMGP/rRjxg/fjzFxcUEAoFuTRxidUx+vz/icwcAFi1a1K/4vV4vRUVF1NfXt0/LzMzE5/NFVavfn331h9fr5frrr+ett97qNS4nr5OwnJwcqqqqKCkpaS/H0tLSmMUWbrbk8XgoLS2lrq6OvLy8TuUw0PflUJ0/6LlcRETcSom+jHher5dt27axbNkyVq5cCQQ/0CsrK2PSzrjrfoqKitoTlqysLEpLS4cleSsvL2fZsmWUlZVRVVWFz+ejuLi4W6ISHie8qKiIH/zgB/h8vvYkr2ttZqyOye/39zhOvM/nw+PxRB1/uIlGpC8O4XHV+xLtvvoj2ricvk4gWGteU1PDypUrycnJae8rEqvYysvLqaysZN26dQAsX76c4uLiTssM5n05FOcPei4XERHXcnpYRje8YjmO/mjXcdx3GbyRVp719fURxzqvr6+3eXl5EcdwH864rrjiiriKa7iFx9HvOF7+YIy06zOe798apzy2VJ6xpfKMHTSOvoiMVOGx3LvWtHq93vYa82g6Gg9VXN/4xjc6TXc6LhERkd4o0ReRuNHbyDDhadnZ2cMaE3wQ1zPPPNNtnpNxiYiI9EZt9EUkbmRlZeHz+cjPz2f58uXk5uYyefJkqqqqWLlyJVlZWRQUFDgW19KlS3n++efjJi4REZHeqEZfROJKZWUlpaWlVFVVkZ+fT3Z2NqWlpaxYscLR5jGVlZX8+te/jru4REREeqIafRGJOwUFBXFZQ37FFVdwzTXXOB2GY/Ly8gj2BRMRkZFANfoiIiIiIi6kRF9ERERExIWU6IuIiIiIuJASfRERERERF1KiLyIiIiLiQkr0RURERERcyDXDaxpjsgAfMAXIAvxAkbU2EOX69UAxsC40yQfk9mcbIiIiIiLxwhWJfijJz7HWlnSYthzYZYzJtNb6o9iMFygNvSD4RSFXSb6IiIiIjERuabrjs9aWdZwQSvrr+CBx70sRkE2wFj/bWhvtFwQRERERkbjjihp9oNAYU2etreoyvQpYHu1GrLV1sQ1LRERERMQZbqnRh2DTGxmhAoEAxhgWLlw4JNuvqKjAGENdXey+y5WUlGCMoarqg++XZWVlMdlPuDxyc3MHG6Z0EOmcuVmsrkcRERmZXFGjb63N7GGWl2DznaiF2vtPjvDrgEgnHo8Hr1ffL0cSnTMRERlNXJHoR2KM8QB5QH6Uq8wxxuQBddbaOmNMMVDfte1/h+0XAAUA06ZNY9OmTVHtZNKkSTQ2NkYZ0ugRLhNr7ZCUT0tLCwB79+6N2faXLFnCkiVLgA/i37dvX8T93HbbbRx77LGcc845UW07vG5ra+ug4h3s+m4T6Zz1R2trKzfffHO/zqWTeroeh0J/r3EYedfnvn37or7XD7empqa4jW0kUnnGlsrTOa5N9AkOlVlhra2IcvnSjrX41toiY0y9McYfqXY/9AWgDGDmzJn27LPPjmonL774ImlpaVGGNHq0trYCYIwZkvJJTk4GYMKECUNa/uPHj4+4nx/84AcsWrSICy64IKrthMsjMTFxUPE2NjbqeouhxsbGfp9LJ/V0PQ6FgZTLSLs+x48fz+mnn+50GBFt2rSJaD+HpG8qz9hSeTrHTW302xljfASH24y2Np8emupUEPzCICIiIiIyorgu0Q812SkCzovB5uoJPnxrRJs+fTrGmPbX9OnTnQ4pahUVFaSnp1NXV0dRURGZmZmkp6eTnx/1d7hOGhoaKCoqIjs7m/T0dHJzcwkEAp2W6akDY3idvpbrKD8/H2MMgUCgfXljDCUlJT2u01+BQKC9bIwxZGZmUlRUFHHZuro6srOzMcaQnZ1NWVkZ+fn5ZGZmtnf8DZe53++nqqqK3NxcCgsL+7Wv/p63aOIK7z8/P5/09PRej7Mnkc5ZtLHm5+czceLEXs9lX/H1VrZVVVVkZ2dTV1dHSUlJe3lEukbD+yosLOwUr9/f94jAsT6Hw3GNi4jIwLgu0QdWAfn9edBVqIlO3tCF5KwdO3b0+nc8a2hoaE+ewglKTk4OFRUV/U7ygPakqbCwEJ/PR1VVFeedF4vvhJEVFxdTWVkJgM/no7KyksrKSgoKCmK2j3Xr1lFVVUVeXh7l5eXk5eVRUlLSnkCG+f1+srOz8Xq91NbWkpOTQ2FhIV6vl9LS0k7lGQgEKC0tJTc3l4aGhvZkO9p99ee8RRuX3+9nxowZ1NXVUVxcTGFhYcR991e0sRYXF3PfffcBkc9ltPH1VLaBQKD9C091dTWLFy8mKyur/QtA1zKbMWMGVVVVFBYWUlxcjN/vJzMzs88RhWJ9DofjGhcRkQGy1rrmRbCZjbfLtKwo1qvtul6H7dX2tf7xxx9vo7V169aol40VoNsr3uzatcsC9pxzzuk0vbS01ALW5/N1mg5Yr9cb9fbLy8stYEtLSztNLygosICtrKzsts/a2tpOy2ZlZVmPx9Prcj2tC9iCgoKo4w2XR9fjjpbP57OA3bNnT/u0goKCTvFba63H47F5eXmdpoXLCrDl5eVR76uj/py3aOMK72fXrl3d9tNxWm96O2fRxLpnz54ez2U08fVWtj1do+HtdrxGfT5ftzKzNniNdoy5p+uxp9g76u97r7/XuLW20/U5Ejhx/47W448/7nQIrqLyjC2VZ+wANbYfubFravRDo+Cstd2fZpsTxeqR1oPgqD3RPllXhlDX8eSzsrJoaGjo93a6Dq1YXBzsglFeXj7w4OJQuLy2bdvWPs3v93c7/pycnB6be+Tl5ZGX1/cPXeF9RdpONOctmrgCgUB7LXRDQwN+v7/TerEYF38w11h/4+utbLuWRWlp8BYUvkbDzX5WrFjRbd1wzX5FRbRjEAQN9hyKiEh8csWoO6HOt5lATWgcfIDJoX+zuyxbDwSstR2nVxljCmyHoTSNMcsBv+1heM2RJCMjo1NznYyMDAejGZisrN67SkRqJx/NmOkejweInOCMJOHkrrKysj3J7Mrr9VJTU9NpWk1NDT6fL+I2Fy9ePOB9hfV13qKNKzy/oqIiYhIbi/MXTaw96W98PZVtJOFrOLyN8LUeKd6cnGC9RnV1da9f0mJ9DkVEJD65pUa/ElhOsAlO+FUZenUVADpVR1lr6wh+SSgOvUpD013xWNLt27d3+hln+/btTofUb5MnT+51/rJly8jOzu70Gkgb/pGooqKCzMxM1q5dS35+PuXl5SxfvrzbcuG21oWFhdTV1VFYWEggEIhYMwzda5b7s6+wvs5bf+OqrKyM+NNkbzFEK5pY+xJtfLF4aFekDrrR1LQPxTkUEZH45IoafWut6cey2T1Mr6OfT9GV+FFbWzug9cI1mcNRazlUzR2WLVvW3rEyLFJTlrq6Onw+HzU1NZSVlZGVlUVlZWW/jj3affVHNHGFa6rDyzqt67kcyvi6XqPhfyPV2odr++fMmdPj9obiHIapSY+ISHxxS42+SFS6NlEI1/pHakrRMWkJj4gyUB6PJ2IN7GAFAgECgUC3GuLq6upuy1ZXV5Obm0ttbS3WWmpra/uVlPZnX/0RTVwejwefz8fKlSu7lWM4ruES6VzGMr6u12h4NJzwNer1esnKyqKsrKzbdouKivB4PD022xmqcwhDd42LiMjAuaJGXyRa4WYh4aEbwx0oO9Yeh5OgcEfd8NCCg5GTk0NVVRUlJSXU19cDH3Sy7I3f76esLHI3kUWLFrX3QygpKSEQCJCbm8vatWsjthP3er0UFRW17x8gMzMTn88XVa1+f/bVH9HGVVpaSnZ2NjNmzGDFihV4PB5qa2spKytrHyZyOPR0LmMVX/ga9Xg8lJaWUldX1+0aLS8v77SvcAx+v799qMtIhuocwsCvcRERGTpK9GVUCCc44TG/V65cyeTJk1m+fHl7Qh/m8/nIy8ujoqKCmpoacnJyKC0tpby8fMCdPouKiqipqWHlypXk5ORE/cXB7/f3OE68z+fD4/FQXl7OsmXLKCsro6qqCp/PR3Fxcbcka8qUKQARvzh0bcrRk2j31R/RxuX1etm2bRtFRUXtSW1WVhalpaXDluRDz+cyVvGVl5dTWVnJunXrACJeo+F9LVu2jJUrVwLBRLuysrLP9v9DcQ5h4Ne4iIgMof6MxanXyBxHfyQZaeNqx7twedbX10cc57y+vt7m5eVFHL99OMRrXD0ZyuszPI5+x/Hy3W6kvd/j+f6tccpjS+UZWyrP2GG0jqMvIj0Ld7bsWsvq9Xrba8wH2qHZjXGJiIi4gRJ9kVGg46gwXYWnZWdHHJBqSMVrXCIiIm6gNvoio0BWVhY+n4/8/HyWL19Obm4ukydPpqqqipUrV5KVlUVBQYHiEhERcRHV6IuMEpWVle0jDeXn55OdnU1paSkrVqxwtHlMvMYlIiIy0qlGX2QUKSgoiMsa8niNazjl5eUR7GclIiISG6rRFxERERFxISX6IiIiIiIupERfRERERMSFlOiLiIiIiLiQEn0RERERERdSoi8iIiIi4kJK9EVEREREXEiJvoiIiIiICynRFxERERFxISX6IiIiIiIupERf4kJdXR3p6ekcffTRpKenR/Xy+/1Ohx1TFRUVZGdnY4whMzPT6XB6VFJSgjGGqqoqp0MRERGRXoxxOgARAI/Hw4oVK9i/fz/jxo1rn15UVNQ+r6vJkycPZ4hDqqKigvz8fAoKClixYkVcf4nxeDx4vV6nwxAREZE+KNGXuOD1elm+fDmNjY2kpaW1Ty8qKmLy5MksX77cweiiV1ZWhtfrxefz9Wu9lStX4vP5KC0tHaLI+q+nYykoKKCgoMChqERERCRaarojEkNFRUWUl5f3e726urq4qyUf6LGIiIhIfFCiLyIiIiLiQkr0ZcQKBAIUFRWRmZnZ3oG1qKio23IVFRXtnXerqqrIzc2lsLCwfX5dXV17J9js7GzKysrIz88nMzOT3NzcTvvLz88nPT29277y8/MxxhAIBCgrK8MYgzGGkpKSXo8h3LEVaF8v3BE3/HddXV2ndbKzs0lPT494jHV1de1lkp6eztKlS3vcd2FhYbeyCx9jb8fSU1yBQKB9m+np6eTn53fra9BTnPn5+b2Wk4iIiPSfEn0ZsdatW0dVVRV5eXmUl5eTl5dHSUlJpyQ+LBAIUFpaSm5uLg0NDe0JvN/vJzs7G6/XS21tLTk5ORQWFuL1eiktLW1P5v1+PzNmzKCuro7i4mIKCws77au4uJjKykoAfD4flZWVVFZW9tmWvaCgoNt6A2ku09DQ0J6kh78A5eTkcO+993b78hMIBMjMzGTdunXt/QLy8vIoKyujpqZmQMcSLp+qqioKCwspLi7G7/eTmZnZaXSenuKsqKiI+CVNREREBk6dcePQ2b+t6zZt0amH8dWPHknzgVY++4dnus2/POdwLp97OO81HSBvzfPd5n/lrA+x+PQM3ty1j8vu3Npt/rc/eTQLTpzKy+/upbDi5W7zv+87Ft/xk3n67Ua+cd8rneZt+mpWfw4vZrp2Cs3Ly6Ouro6ysrKInVpLSkravxCEFRcX4/F42pPr0tJS1q1bh9/v79QJtbCwkEAgwLZt2/B4PEBw9JlwUuv1etvb2PenM67H42lftuP/Byr8BQWC5WOMoaKiguLi4vZlioqK8Pv91NfXd+oX0HGZ/h5L+AtPfX19+7SCggKys7MpLCzsND3aOEVERGRwVKMvrtKxpr6rvLy8Tkl+eLmunWBzcnI6rR8IBNp/OWhoaMDv93daL57Gk+/Y1AjgtNNOo6GhodO0srIy8vLyYtb5N9wkKtIQqOGa/YqKil7jzMrK6haniIiIDI5q9ONQbzXkKWMTe50/NXVsr/OPSh/f6/yZh03odf5pH0pzrAY/knASWVlZ2Z6A92Tx4sXdpnm9XmpqajpNq6mp6VSLHZ5fUVHRLWENxxAvsrJ6PzfhdvVdE+3BCG8z0r5zcnIAqK6u7vQlq684RUREZPBUoy8jVkVFBZmZmaxdu5b8/HzKy8t7HW8/Ug12uK14YWEhdXV17U10ItVOV1ZWYq3t9oqnMf6dfIhYIBDoNq2nWno3PexMREQkXinRlxFr2bJl5OXlUVtbS0FBAVlZWUyZMqVf26irq8Pn81FTU0N2djY1NTVUVlZ2qnEO10p3HWWmJ/HcBCV8XLW1tVEtH82xhLdZXV3dbV64zObMmRNtiCIiIhIjSvRlRAoEAgQCgW619JGSzd5UV1eTm5tLbW0t1lpqa2u7dT4Nd5JduXJlt1rrcBwdl41Usz0YHZPtQCAQ9ReOnoRH2Ona5CgQCHSaFu2xeL1esrKyKCsr67Z8UVERHo+nW98IERERGXpqoy8jksfjwev1UlJSQiAQIDc3l7Vr10ZsQ98br9dLUVFRp1FhMjMz8fl8nWr1S0tLyc7OZsaMGaxYsQKPx0NtbS1lZWWdRvLJycmhqqqKkpKS9m1GGgEo2tjgg9FwwsNRDtaqVavanx1QUFBAZmZm+7EUFxe3N0Xqz7GUl5d3Kp/wsn6/v32oThERERleSvRlxCovL2fZsmWUlZVRVVWFz+ejuLi4X4l1uKlPWVlZt3nh8fkhmHRv27aNoqKi9gQ2KyurfQz6sKKiImpqali5ciU5OTmDSsx9Ph95eXlUVFRQU1NDTk4OpaWllJeXD6oDcPhLSlFRERUVFe3H0jHJ7++xhMtn2bJlrFy5Egh+UaisrIzZ6D4iIiLSP8Za63QMI97MmTPtyy93H3s+khdffJFZs2YNcUQjV2NjI2lpacOyr/ADnQoKCjp9OfD7/e1JcGlpaZ8PvYpnw1meo4HKM7ZGWnnG8/1706ZNnH322U6H4Roqz9hSecaOMabWWpsT7fJqoy+jVnj8+6411V6vt70mP9pOqyIiIiLxRom+jFq9jaYTnpadnT2sMYmIiIjEitroy6iVlZWFz+cjPz+f5cuXk5uby+TJk6mqqmLlypVkZWWN6GY7IiIiMrq5JtE3xmQBPmAKkAX4gSJrbSDK9T3ACiA8/EqmtXbwQ5xIXKusrKSsrIzS0tL24SG9Xi8rVqyIqwdhiYiIiPSXKxL9UJKfY60t6TBtObDLGJNprY1miJJyoDC8rDHGa4yptNbmDk3UEi8KCgpUcy8iIiKu45Y2+j5rbafxEUNJfx3Q51iLxpg8wN/xC0GHhF9P+hERERGREcctiX6hMcYXYXoVweY8fVkMRBpepRIoHExgIiIiIiJOcEuiDzCYp/L4CLbp78oPRD1WqYiIiIhIvHBFG31rbWYPs7wEm+/0KNQJ1wM0RJgdCM0TERERERlRXJHoRxJK4POA/D4WnRzNtqIdvUdEREREJB64NtEHioEKa23FUGzcGFMAFABMmzaNTZs2RbXepEmTaGxsHIqQXKG1tVXlE0Mqz9hSecbWSCvPffv2RX2vH25NTU1xG9tIpPKMLZWnc1yZ6Ic65uZYa4fssaahUX7KAGbOnGnPPvvsqNZ78cUXSUtLG6qwRrzGxkaVTwypPGNL5RlbI608x48fz+mnn+50GBFt2rSJaD+HpG8qz9hSeTrHTZ1xgfYmO0XAeVGuEqltfidqtiMiIiIiI43rEn1gFZAfbXIeWi5A5FF7vKF5IiIiIiIjiqsSfWNMMVDUMckPPTW3LzVE7pSbSXAsfhERERGREcU1iX6oc+zajk+3DYlmHPxyIDfCdB+wdrCxiYiIiIgMN1ck+qHOt5mh/2eFXr7Q9Owuy9YbYzo9BTfUsdZrjPF2WC4LaBiqUXuks0AggDGGhQsXDngbJSUlGGOoqoruR5j+Lt9RWVkZxhjq6np9TEOvcnNzyczs6REQIiIiIoPjikQfqASWA7UdXpWhV1cBInfAPQ8oNMYUhH4dWGytjVTLP2JYa9m8eTP5S77AhIkeEhITmTDRw6JLLmPLli1Ya50OMaY8Hg9eb/QPSO7v8iIiIiIjiSuG17TWmn4sG3HIzVC7/qJYxeS0gwcPsvTKZdy/sYp9s+fTtqQMUtJpbt7F3S89yob5n2fBPB9rVq8iKSnJ6XBjoqCggIKCgm7Ty8rK8Hq9+Hy+qJYfKXo6LhERERFwT42+dGCtDSb5m1+ieUkpbTmLIXUqJCRC6lTachazd0kp921+kaVXLnNdzX5XRUVFlJeXOx1GzLn1uERERCQ2lOi70JYtW1i/sYrmeddBUnLkhZKSaZl3Pes3VlFdXT28AYqIiIjIkFOi70I3/foWWmbP7znJD0tKpmX2fG66+ZbhCWwAKioqSE9Pp66ujqKiIjIzM0lPTyc/P7/bsl07yObn52OMIRAItM8zxlBSUhJxeQh2Cg7vxxhDZmYmRUWDb9FVVVVFdnY26enpZGdnU1FRQUND964i0ey/r+PquI2JEyfG7BhERERkZFGi70IPPrCethOiezBw28xzeXD9+iGOaOAaGhoIBALk5+e3J7A5OTlUVFT0mbwWFxdTWRnsj+3z+aisrKSysrLXdvnr1q2jqqqKvLw8ysvLycvLo6SkhMLCwgEfQ11dHbm5udTV1VFQUEBhYSErV66MOGJPNPvv67g6bmPNmjUxOQYREREZeVzRGVc6a9nbCCnp0S2ckk5LU+PQBhQDXq+X0tJSINiJ1hhDRUUFxcXFva4THlUn2k6rXTvo5uXlUVdXR1lZWfv++2vZsmUA1NfXt8dTUFAQcWjNaPbf13F13EZjYyOXXXbZoI9BRERERh7V6LtQ8oQ0aN4V3cLNu0hOTRvagGIgN7fzSKdZWVkRm74M5b79/q7PYutbIBBor8nvOpSnx+MZ8v3HchsiIiIysqhG34Xmn7+Au196NDjaTh8SXn6M+QsWDENUg5OVlTVs+/L7/VRUVFBZWYnf7x9Uchxetz8PxorF/sPb2LhxI6+//roSfBERkVFINfou9O1rryZ564NwsKX3BQ+0MH7rA3z7mquHJ7BBmDx58rDsp6KigszMTNauXUt+fj7l5eUsX7580NuNtvY+FvvvuI0LL7wwZscgIiIiI4tq9F1o7ty5LJjn476NN9Ay7/rIo+8caCH54Ru4YF4uc+bMGf4gh1m0zXyWLVvW3hE2rKqqasD7DTfXqa2tHZL9RzqujttobGwkLS1tUMcgIiIiI5Nq9F3IGMOa1atYeMYsJtxZSELNWmjcCa2HoHEnCTVrSbmrgIVnzGLN6lUYE/WDhUckj8dDIBDoc7lAIEAgEOjWln4wzxnweDxkZWVRVlbWKYZAINCtOU1/9x/puIbiGERERGRkUo2+SyUlJXHHmtuorq7mF7/6DRvWXkVLUyPJqWnMX7CA7/zi3lFRkw+Qk5NDVVUVJSUl1NfXA0Qcfcbj8eD1eikpKSEQCJCbm8vatWupqKgY1P6Li4vJzc1lxowZ7aMEFRcXEwgEOjVJ6u/+ezqujtv42Mc+xv333z/oYxAREZGRRzX6LmaMYe7cuay748807d5Fa+shmnbvYu1f1oyaJB+gqKgIj8fDypUr8fv9ER+2FVZeXt5eA19UVMTkyZMpLi7uVkPeH+Gx7r1eL0VFRVRWVlJaWkpeXl637fZn/z0dV8dtXH/99TE5BhERERl5jLXW6RhGvJkzZ9qXX345qmVffPFFZs2aNcQRjVzhNuUSGyrP2FJ5xtZIK894vn9v2rSJs88+2+kwXEPlGVsqz9gxxtRaa3OiXV41+iIiIiIiLqREX0RERETEhZToi4iIiIi4kBJ9EREREREXUqIvIiIiIuJCSvRFRERERFxIib6IiIiIiAsp0XeAnl0gIjKy6L4tIiOREv1hlpCQQFtbm9NhiIhIP7S2tpKYmOh0GCIi/aJEf5iNHz+e5uZmp8MQEZF+aGpqIiUlxekwRET6RYn+MEtNTSUQCOhnYBGREaK1tZWGhgYmTpzodCgiIv2iRH+Ypaenc+jQId555x3279+vhF9EJA5Zazl06BCBQIDXX3+dCRMmkJaW5nRYIiL9MsbpAEabhIQEjjrqKBoaGnjjjTc4dOiQ0yHFlX379jF+/Hinw3ANlWdsqTxjK97LMzExkZSUFKZOnUpaWhrGGKdDEhHpFyX6DhgzZgyHHXYYhx12mNOhxJ1NmzZx+umnOx2Ga6g8Y0vlGVsqTxGRoaWmOyIiIiIiLqREX0RERETEhZToi4iIiIi4kBJ9EREREREXUqIvIiIiIuJCSvRFRERERFxIib6IiIiIiAsp0RcRERERcSEl+iIiIiIiLqREX0RERETEhZToi4iIiIi4kBJ9EREREREXclWib4zJM8aUD2C9emNMgTHGE3rlGWNKjTGeIQhTRERERGTIjXE6gFgwxpSG/usFJg9gE16gNPQC8AO51trA4KMTERERERl+rkj0rbWFAMaYAqBwAJsoAqoIfklosNbWxTA8EREREZFh54pEPxaU3IuIiIiIm7iqjb6IiIiIiAQp0e/AGJNljPE5HYeIiIiIyGAZa63TMcRMuI2+tTa7n+uVA2uBOmut3xhTDNRba8v62FcBwLRp07LXrVs3iMglrKmpidTUVKfDcA2VZ2ypPGNL5Rk7KsvYUnnGlsozds4555xaa21OtMurjX5QqbW2KvyHtbYoNOSmv+P0jkJfAsoAZs6cac8+++zhidTlNm3ahMoydlSesaXyjC2VZ+yoLGNL5RlbKk/nqOkO0EMyXwEUD3csIiIiIiKxoES/Z/VAltNBiIiIiIgMxKhP9ENNdPKcjkNEREREJJZGfaIPBIBIY+hn9jBdRERERCTuKdGHtdZaf4TpeUDpcAcjIiIiIhILbkv0PaFXRKFmOrVdJleFhsrsuNxywN/b8JoiIiIiIvHMFcNrhsa99wCLAE9oXPwGgsNmdmx+EwhNb2etrTPGhLdBaDv11trcoY5bRERERGSouCLRt9YWhf5b2MdyER+kFfoyoPb4IiIiIuIabmu6IyIiIiIiKNEXEREREXElJfoiIiIiIi6kRF9ERERExIWU6IuIiIiIuJASfRERERERF1KiLyIiIiLiQkr0RURERERcSIm+iIiIiIgLKdEXEREREXEhJfoiIiIiIi6kRF9ERERExIWU6IuIiIiIuJASfRERERERF1KiLyIiIiLiQkr0RURERERcSIm+iIiIiIgLKdEXEREREXEhJfoiIiIiIi6kRF9ERERExIWU6IuIiIiIuJASfRERERERF1KiLyIiIiLiQkr0RURERERcSIm+iIiIiIgLKdEXEREREXEhJfoiIiIiIi6kRF9ERERExIWU6IuIiIiIuJASfRERERERF1KiLyIiIiLiQkr0RURERERcSIm+iIiIiIgLKdEXEREREXEhJfoiIiIiIi6kRF9ERERExIWU6IuIiIiIuJASfRERERERF3JVom+MyTPGlA9gPY8xptgYUxB6FQ9FfCIio5m1ls2bN5O/5AtMmOjh3HPPY8JED4suuYwtW7ZgrXU6RBERVxnT3xWMMZ8HAtbax4YgngExxpSG/usFJg9gE+VAobXWH9qe1xhTaa3NjVWMIiKj2cGDB1l65TLu31jFvtnzaVtSBinpNDfv4u6XHmXD/M+zYJ6PNatXkZSU5HS4IiKuMJAa/auAKmNMqzGm2hiz0hhzrjFmYqyDi5a1ttBaW0gwYe8XY0we4A8n+aHt+TvMExGRQbDWBpP8zS/RvKSUtpzFkDoVEhIhdSptOYvZu6SU+za/yNIrl6lmX0QkRvqd6FtrPwVkAouBWiAfqAJ2GWNeMcb8zhjzudiGOaTCx9FVJVA4zLGIiLjOli1bWL+xiuZ510FScuSFkpJpmXc96zdWUV1dPbwBioi41IDa6Ftrt1lrK6y1V1lrjwPSCSbMjxFMju8O1fj/1sma/ij5AH+E6X4gZ5hjERFxnZt+fQsts+f3nOSHJSXTMns+N918y/AEJiLicjHpjGut3R1K/AuBbILJ/iLgOGCbMebUWOwn1owxHsADNESYHQjNExGRfrLWcthhh2GMoXzdOtpOOC+q9dpmnsu6u+5i+vTpQxyhiIj79bszbl+stf82xiyy1q4iWLPvAyqMMbnW2tdivb9B6rPjrjHGY60NRJheABQATJs2jU2bNsU8uNGoqalJZRlDKs/YUnkGff7zn2fXrl3tf0+aNIlvfHcFz7wV4NWGA/x3fxJ7dr7DoZ07gwu0HYKU9Og2npIOba3s8H6azMxMjjnmGJK92Zz4IQ+zZ3yID33oQyxevLjT/tPT0/nrX/8ay0MccXRtxpbKM7ZUns4ZyKg7NwLLAQtUAHdZa+/paXlrbZUxZg5QBKwYaKDxxlpbBpQBzJw505599tnOBuQSmzZtQmUZOyrP2FJ5QmNjI7v2WZhxJmx7CoDdn/w2P2qYAykEX60HgQ7t7BOToHlXsANuX5p3QdJ42LsLv9+P/8134PTvs6EJ2LwLGp6H0y+D/2yCt54FYNeuXaP+vOjajC2VZ2ypPJ0z0Br9TCAXyCNYa2+BOoLt2j0Em720s9YGjDGR2sGLiEic27lzJ7dVbOD2f73KiwenwlXlcHA/3HJ+cIFtT3Fc6iFmThvP3MwMPn7KcZx84lVMu++6DzaytRLmLul7Z1srg78APLueLVu28MzzW3l4ay1btzfxZksCjUnpcMJ58P5rwUR/wmQ4eT7fKf4tX82fh9frHZIyEBEZiQaU6FtrtxGszS4zxkwi2BHXR7BTbp219nvQXvs/0Vr7VYJfDuJNpLb5nURqtiMi4nZvvPEG99xzL/fc81f+1jgFe/bX4LCjYKcfnvozvPrPDxZ+bgOvPPtgt21kZGSwY8cOOLQfasvh9At775B7oAXqKuDQATIyMpgzZw5z5szhyx0WaWlpISUlBUyoi9kRJ8JZS7np3QRu+uk/OWz3bVwwezJfz/dx8kknYYyJSXmIiIxEAxle83vGmO8aY74c+nu3tbbMWrvIWvupcJIf4gcuNsb8Dlgbo5hjJpTEBwg+aKsrL11+mRARcavp06djjMFMPwHz8S9zzA8e4xu//DNPPPEEiW8/y8wdm/jxjNfZ/vMFZPgfhp317etmZGRE3Ob27dux1tLW1sbFF11I8sYb4GBL5AAOtJD88A1cfNGFtLW1sX379oiLJScnB/dn24ITXvk7aXd8iZxdfyNxz3bePfJj/GHvqZx61jkkJSVhUqdgEsdgjFEHXxEZdQZao/9zY8wMY8xp1tqne1muvR17HKshcqfcTILPBxARcbV/PrmZHYdlw2fyYfLR0NYKbz7NJz56FoU/+ybz589n0qRJ7cv3lIT3xBjDmtWrWHrlMtbfWUjL7Pm0zTw32PG2eRcJLz/G+K0PcMG8XNasXtVnLXxP+9+3bx/3P/wof3zoH9QlG97b2wq534aM4+Hf97Dj2fW0tbWRkBCTAedEROLegEfdCTXfcYNygv0Nun4h8QErhz8cEZGh19bWxgMbHuKmn5fwt3/8C770l2BH2IdLoP5J2LeHJ2L4hNqkpCTuWHMb1dXV/OJXv2HD2qtobtpDSupE5i9YwHd+cS9z5swZ1D7Gjx/PooXzWbRwPodu+QFJSUnwzP1w2oXwsS/B3EuYmv8Dvpd7HNdecTHjxo2LzcGJiMSpmA+v6TAPvYx9b4ypBwLW2uzwNGttmTGm0Bjjtdb6Q8tlAQ3W2oohjldEZFgdOHCA3/xpLT97+GUa0jLhH/9iUtoEdt95NTS9N6T7NsYwd+5c1t3xZ2BoR+IYMyb08bZtc/A11Qs5+eyaeS5Fq//Er25YwTXXXMtVVxXi8XiGJAYREae54vdLY0yxMaaU4PCdXmNMuTGmNJSwdxQgcgfc84BCY0xBaHz8xdba3KGNWkRk+OzevZtv/ewWPEtu5DtbM2g4+pOkNL3Nj278BW+88QYZExI7Ld9Tu/uRpNMxvOcn499/4lenvMfsQ6/yzjvvsOJPDzO1cDWLvnMjb775pnOBiogMEVfU6Ftri0L/LexjuewepgcIjvMvIuIa06dPD456A5AxEy79LRzRwpS3/sH/fuYEvla8nLFjxwL9b3c/EvR0TNdcsYRHHnmEb67awIueoyi3Uyj/nw1QvRb+8wTYNjIyMlxZJiIyuriiRl9ERDo7cOAAOw4lw8xzghN2vAxVv2JtbgI7y3/EN790SXuSP9oYY/j0pz/N1opf868vzWDOridgzDiY/31Y+GOAD74giYiMYK6o0RcRkQ/cde+DXLVmM1xWBnvfh1f/EXxa7bPrWXTB/U6HF1fOmpvDlrk5GJMAH/44HNwXnDFmLD/5zWpWfPWLJCYm9r4REZE4pRp9EZE4Zq1l8+bN5C/5AhMmekhITGTCRA+LLrmMLVu2YDuMjPP81hc58dL/YckjB9l9zMfh2fXw58Jgki99sPDK3+C1LcE/T7+I6149gsMXXc9Dj27qvGQ/zomIiJOU6IuIxKmDBw9yydIrOPf8i/jrfyfQvKQMe81GmpeUcffbKZw7//NcsvQKdu7cybe+9S1OO+9CtmacS2LDa3zn8FfJeGEd7NvTvj03dLAdKl3LJu2tJ5nwVjU7jz2Pz67bydzLv8+2116L+pwcPKgvVyLiPDXdERGJQ9Zall65jPs3v0TzklJISv5gZupU2nIWs/fUC7h7/fWUn/JRWre/gjGGzzVX8rtff4uMjAx+vvwrzh3ACBOp421zczPfuPH3/GHHeKozzuW4gpuZ9e4T+JvH0NLLOblv4w0svXIZd6y5rc+Hf4mIDCXV6IuIxKEtW7awfmMVzfOu65xQdpSUzMEFN9DavIcTs8+gtraWv/6+WDX3MZKSkkLZDd9i28/O58zAJtpqKnjhP/W0zLu+13PSMu961m+sorq6engDFhHpQom+iEgcuunXt9Aye37PCWVYUjKcfD4nzjye008/fXiCG2WOOfponvzjDZyTcyJk50d1Tlpmz+emm28ZngBFRHqgRF9EJA49+MB62k44L7qFZ+fy4APrhzYgYfNTT8Ls6J6l2DbzXB5cr3MiIs5Soi8iEoda9jZCSnp0C6ek09LUOLQBic6JiIw4SvRFROJQ8oQ0aN4V3cLNu0hOTRvagETnRERGHCX6IiJxaN5nPgsvVkW1bMLLjzF/wYIhjkjmn7+AhJcejWrZhJcf1TkREccp0RcRiTObN2+m6r/AM+vhYEvvCx9oYfzWB/j2NVcPS2yj2bevvZrkrQ9GdU7a6u7FO3P28AQmItIDJfoiInFg+vTpGGMwxnDmmWey5+1tJE0+nHEP/rDnxPJAC8kP38AF83KZM2fOsMY7Gs2dO5cF83wkb7yh13Ni7r8Ojj6NG9+ZiZnlaz+v06dPH96ARWTUU6IvIhIHdux8Hz5eAHMuDk7Y9hTvP/sEn/vIiUy4s5CEmrXQuBNaD0HjThJq1pJyVwELz5jFmtWr9GCmYWCMYc3qVSw8Y1av52TxJ07mf6/8HOx6Cz72JRgzDoAdO3Y4fAQiMtroybgiIg77xaq/wOJfweGzoO6v7dPT0tK4Y81tVFdX84tf/YYNa6+ipamR5NQ05i9YwHd+ca9q8odZUlJS1OfkpwmJkHYYHNoPCWNg6gwOHTrEmDH66BWR4aG7jYiIQ5qamjj/m8U8MXYOpE+E9T+CV/7WaRljDHPnzmXdHX92KErpKupzYttgz/bg/+dcDGct5cOXXMcTN32Vo486augDFZFRT013REQc8Mwzz3DKJ+bxROonSNj9Nqn3frtTkp+RkeFgdBILnc7hv++BV//Ba0flctx3y/l/9zzoXGAiMmoo0RcRGUbWWn5+SylnnHEG2/79T46t+z21RZ+g8e1XsNa2v7Zv3+50qDJI27dv/+Cc7m9i+x++yuz/PsLBjNl8YWMzF3/7Jxw4cMDpMEXExZToi4gMk0AgwNzLv8/ylz/E/sNPoaCggK0b/8JpJ5/odGgyDDIyMnju//2Mb2a8Cs0NrL39j3zsYx/D7/c7HZqIuJQSfRGRYbDpn5s58su3UDP1PBLf8/PbG75LaWkpycnJTocmwyghIYH/+5+r+ddXT+GYVEt1dTWzl/6Q3//lbqdDExEXUqIvIjKEpk+fjjnsOM75w0vsPfpMzOY/s/WGz/KVy/KdDk0cdNZZZ/Hvf/+bT+UvZf/pi/nKvxIw3rM05r6IxJQSfRGRIRIIBIJjpx8+C8YmQ/l3sP/8E8d/+DinQ5M4kJ6ezsa1f+I673+D4/F/7qfBZymYBI25LyIxoeE1RUSGwObap7noqu8G/3j2AXj5cdi/19mgJO4YY7jhm8v48Zix8MmvwJzFkH4k3H+906GJiAuoRl9EJMZ+ser/cdavt/D2GV+H8WnBiUrypTetB+Gxm+GhlfDCwwBcd911tLa2OhyYiIxkSvRFRGJk//79fPZrP+K7T6di0zLwHajmsImdO9tqfHyJpP26eLEK6v8JwE/WP8eJX7iOnTt3OhiZiIxkSvRFRGLg9ddfx/uFG3ho/Mcxje/ysxN3UVn2Y3bs2KHx8aVPncbct5bKyirGHjeXl4/wccxX/8Cmf252OkQRGYGU6IuIDFJlZSXZ2dn8d1czE177O3/76mms+OrlToclI5jPdx6v/PIyjnr9YVqOzOHcVVv54c1/xFrrdGgiMoIo0RcRGaC2tja+8qNf8alLr+L999/n00kv81rp1/jYmXOcDk1c4OijjuLVO37MBYf+hU2eyI/+M42LryikubnZ6dBEZIRQoi8i0k/Tp0/HGEPiaQv5fcMJ8Mmv8IMf/IAHH1jP1KlTnQ5PXGTs2LHc95sf8Ou5B0j65yrW3b6KCRMmaLx9EYmKhtcUEemnHe/vgk99F06aB69tgQ0r+WH5bqfDEhe75oolnJNzEqec8gAcdTqccSls+KnG2xeRXqlGX0SkH2697Q5YfHMwyX9yDdzzv7Bvj9NhyShw8sknB/8zPi34ELYv/B4On60hOEWkR0r0RUSicODAAW6++Wau/vJSaNoZTPCfvB1sm9OhyWjzyt/gzq/DoQOw6P84Zel17N6tL5si0p0SfRGRPrz93/9ywiX/yz0PPUpSYgITH/85bHuqfb7Gxpfh0n6tveeH//cVeL2OrdN9XPHj3/P00087GpuIxB8l+iIivXh40z/I/MZf2HbMZ0g54yL+9re/sXv3bo2NL47oNN7+vkb8v1zCsS/8hV21D3HWWWdx+5//4nSIIhJHlOiLiERgreV/bipj3pp69h9xGplvbmTN18/nzDPPdDo0kXYzjj2Wrffcyrx589g3Lp3LN1kuuObHHDhwwOnQRCQOKNEXEemiubmZT31pOStf+xCMncCixC28dMfPmDJlstOhiXSTnJzM8uXL+fGPfwyH9rN+zFl8+NIf8tbbbzsdmog4TIm+iAgfjI1vjGHChAlU3bWKxDdr+O1HLWtv+j5jxmg0Yolfxhi+//Uv8fiyE0l+4yneONrHUV+7DTMuVWPui4xiSvRFRCA4HnmyBz5RCAmJ0LKbp39yEV+5LN/p0ESidvbHzsJ/y5fJfGMjHDs3ON5+iMbcFxl9XFNFZYzxACuA+tCkTGttUT/WrweKgXWhST4gFyiy1gZiF6mIxJvW1lY4fDYs+EFwjPL/PAHbX+Kkk05yOjSRfps+fTov3fkzko46Gd7bFpyYlAwHW5wNTESGnWsSfaAcKLTW+gGMMV5jTKW1NjfK9b1AaegF4AdyleSLuNu7777LJ679FSz6JTTugDv/B3bW97meSDwbM2YMbH8p+EfSeLjkVnjzaZ5+7gVOO/lEZ4MTkWHjiqY7xpg8wB9O8gE6JPx5UW6mCMgmWIufba3N7Lg9ERlZrLVs3ryZ/CVfYMJEDwmJiUyY6GHRJZexZcsWrLU8+eSTZF5RwstH+GDb5uC45KEkX2Pjy0jXfg23Hgw+9+G0hWQX/41bbg/+cB3Ne0RERja31OgvBiojTK8ECoGKaDZira2LZVAi4oyDBw+y9Mpl3L+xin2z59O2pAxS0mlu3sXdLz3KhvmfZ6b3KJ6praY1/RiOSU7gH7dew5FHXud06CIx0/H5Dk1NTXzmmpX8w3MmX998iIe3/JCUQD0PPPJ4j++RBfN8rFm9iqSkJAePQkQGwxU1+gTb00eqffcDOcMci4g4yFobTPI3v0TzklLachZD6tRgB9vUqbTlLGbvklLq3tlHK4lce+kCXrnzpxx55JFOhy4yZFJTU/nbH3/C9cdth6b3eWDjI5T//ble3yP3bX6RpVcuU82+yAg24hP9UCdcD9AQYXYgNK8/28syxvgGG5eIOGPLli2s31hF87zrgh0QI0lKhgt/wri0dC655BLVWMqoYIzhR99Yxm8/bjA7XsJe+NNe3yMt865n/cYqqqurhzdQEYmZEZ/oA30+wSb0ZaAvc0Lt+QPW2ipjTLExpmDQ0YnIsLrp17fQMnt+zwlMWFIyB09ZyE033zI8gYnEiccf34TJzo/qPdIye77eIyIjmBnpP8kZY7wEh9TM7trGPlQzXwmk9zV6jjHGZ62t6jKtnuBIPlURli8ACgCmTZuWvW7duq6LyAA0NTWRmprqdBiuMRrLc975F7D/0j8EmyL0pXEn4+8o4KEH7otq26OxPIeSyjN2+lOWQ/kecQtdm7Gl8oydc845p9ZaG3WzdLd0xh20SMk8wU68xQRH4+m6fBlQBjBz5kx79tlnD2l8o8WmTZtQWcbOaCzPAy17ISU9uoVT0jnQsjfqMhqN5TmUVJ6x05+yHMr3iFvo2owtladz3NB0J1Lb/E4GMRZ+PZA1wHVFxAHJE9KgeVd0CzfvIjk1bWgDEokzeo+IjB4jPtEPJfEBgg+86sobmtcrY0x9P8bbF5E49ul5n4EXI/1A113Cy48xf8GCIY5IJL7MP38BCS89GtWyCS8/qveIyAg24hP9kBoid8rNBKL5xA8AkcbQz+xhuojEoS21/6by3WR4Zj0cbOl94QMtjN/6AN++5urhCU4kTnz72qtJ3vpgVO+Rtrp7OfHk04YlLhGJPbck+uUEn2jblQ9YG8X6a3t4Cm4eUDqYwERk6FlrueWWW/hY/jKaTl/C2MmHM37Dj3pOZA60kPzwDVwwL5c5c+YMb7AiDps7dy4L5vlI3nhDr+8Rc/91cPRp/KD+Q+R/8yfs379/eAMVkUFzRaIf6hjrDY3AAwTHwwcarLWdnoobaqZT22UTVV2H0jTGLAf8oW2LSJxqaGjgvMVf4utf/zoHt9WyaN8jvPv041x41mwm3FlIQs1aaNwJrYegcScJNWtJuauAhWfMYs3qVRhjnD4EkWFljGHN6lUsPGNWr++RxZ84ma9fPB8O7aOCM/hw3nd45ZVXnA5fRPrBTaPunAesCA2JCZBprY1Uyx+gSwdea22dMQZjTHFokgeo72F9EXHY9OnT2bFjBySNh3OvgVlLSPU+zeobV5Cfnw/AHWtuo7q6ml/86jdsWHsVLU2NJKemMX/BAr7zi3tVky+jWlJSUtTvkc/94yku+OXDvFlTyfHHdx5TPyMjg+3btztxCCISBdck+qFOuUVRLNdtqMzQ9DrUHl9kRNixYwdMy4T534f0I+Gpv1D3SAUfzvygT74xhrlz57Lujj87GKlI/Ir2PXLOx87krZNnUVCwlXXrXoaPXA7vvATbngq+F0Ukbrmi6Y6IjB4vvfQSZOfDJbfC2BQo/w48eXunJF9EYmvSpEncddddMGYszDgTPvdT8H2z76frioijlOiLyIjQ1tbGr371K04//XRISIRX/wlrlsFbzzgdmsioYIyBQwfgrq/Dljvh5M/CZWXkXXsDzc3NTocnIhEo0ReRuOffto2Zi4v45m/WsW/fPsY9dy88+GPYtwcIthMWkaGXkZEBrQfhH3+Atd8ALHe3ZXHK3I/y1FNPOR2eiHShRF9E4pa1luLf3cbx37uXV4/+DONP/hT33nsv+1pasNa2v9QZUGR4bN++/YP33tvP87dlszjm32XUv/A0H/noR7lqxU84cOCA02GKSIgSfRGJS2+//TanXbKc770whdbDPkzW+5t4Y9XXWLhwodOhiUjIx8/M4aWqdXz3u9/FzjyH0n1ncNTi66j599NOhyYiKNEXkThjreWuu+7ihE8t4dkjPkNiw+v8MmsvNbf9iGnTpjkdnoh0MX78eEpKSthwyw9IfWsL73o/zZz/28y1N/ySQ4cOOR2eyKimRF9E4sL06dMxxpDgOZwlS5bQtPXvnP7mPWwr+RzfuGKxHmwlEuc+c+7Heee2a/n0gSchdRo3vz+LpKwLCT2nhunTpzsdosioo0RfRBzX1tbGjr2tcP718MXVMOlwAGrX3sxRRx7pcHQiEq3U1FQ2/uZ/ueuzqfDaFtj7wfMpNea+yPBzzQOzRGRkeva55/ncj24PJviJY+Cpv0DjuwCqxRcZoRYv+BQXX/DpDybMWQLJE1m/sYoF83zOBSYyyqhGX0Qc0dzczPLv/Q+n3vh3/Md8Bt55AW7/Emy5A9panQ5PRGJpwmTIWcQFd+/kk1euUO2+yDBRoi8iw+6+BzZw4okn8vPilfDaZnz7n2TaP34Ju99pX0Zj44uMbJ3ew5tuJfm+IszBFv6WnsuRX1vNz37zB9ra2pwLUGQUUKIvIsPmrbfe4qylRVx4fyOv7U/h1FNP5cmbrqLylv/l3R07NDa+iIt0GnPfWppfrebF//0kM3c8waHpJ/G/PynhIx/5CM88o6dbiwwVJfoiMqTCo+kYzxEcdc2feWrap0nYt4dvXXM1NTU1nHnmmU6HKCLDZOaHM3nxzz/kTx87wBFj9rJ582ZOKyjBHDFbo/OIDAF1xhWRIbVjxw44/fPw8S8H294/fiuv3vNLZhxztNOhiYgDjDF88eKL+Nxnc1l+3Y8obc6Bj14Bzz4A/7xN7fdFYkg1+iIyJF7YupUFF1wQ/GPMWNi2Gf50Bfz7r0ryRYSJEyfy+1/fFOyEX/dXOHk+fOnPMOdiHnlsk9PhibiCavRFJKa2b9/Ol2/4LQ82Hgn/aQpOrF4LWEfjEpE4dbAFnvgdPLch+MvfGV/g0ws+z4LzPsaNN97I7NmznY5QZMRSoi8iMbF3716+W/I7yl6E1qPOhoPbmec7m9qG59m5c2f7chpNR0Q6ysjICDbXaXgd7ruOlIxjSEg4yPr161m/73g+e1wKf7z+q2q7LzIAarojIoNy6NAhVq1aRcbnivhd42m0TjuOE9/7G09fexoP3XI97777rkbTEZEedR2dZ+/213j11VdZ+pVvwtGnsyH5bI685i985fpf0NTU5HS4IiOKEn0R6cRay+bNm8lf8gUmTPSQkJjIhIkeFl1yGVu2bMFa+8FIOuNSSUqeQEFBAXu3PcP0t//O+gs9PH/7Dzn1JP3cLiIDk5GRwe2//T+e+WYWJ733N1qnHcfvm04jLf8nmOSJnUboieaeJTJaqemOiLQ7ePAgS69cxv0bq9g3ez5tS8ogJZ3m5l3c/dKjbJj/eRbMO48d7wcgOx/mXAy1FVB9J3f++GssWrSIhATVH4hIbJxy4iyeu/2HPFj1BF9e9RjbD58Nh/YHZ5oEduzYwSVLr+jjnuVjzepVJCUlOXswIg5Qoi8iQLAmf+mVy7h/80s0LymFpOQPZqZOpS1nMXtPvYC77rsOjvsYfKIQ3qiD12sAuPjiix2KXETcbr7vk/z3vE+QMGYstB2ChDFw2SqovIm7//k8B3u5Z9238QaWXrmMO9bchjHGuYMQcYCq3kQEgC1btrB+YxXN867r/IHZUVIyLPwxvP0srFkGdy+Hd18Z3kBFZFQyxgSTfICk8bD1YWh6j4Pn/6jXe1bLvOtZv7GK6urq4QtWJE4o0RcRAG769S20zJ7f8wdmWFIynHoB7H67fZJG0hGR4dB+r9nfBP++B045P6p7Vsvs+dx08y1DH6BInFGiLyIAlJevo+2E86JbeFYuE1JSNJKOiAyrjiP0pCSPh9m5Ua3XNvNc1q1dO8TRicQfJfoiEnToIKSkR7dsSjotTY1DG4+ISC9a9jb2657FoYNDG5BIHFKiLzKKtbW1Uf7AI2R/4XuQOAaad0W3YvMuklPThjY4EZFeJE9I69c9i6Rx3HHHHezfv39oAxOJI0r0RUah9957j6/89PdMWlLMoso26jI+BWNTYGtlVOsnvPwY8xcsGOIoRUR6Nv/8BSS89Gh0C79YBWNTufTSSzn8pDP53PL/4z+v+oc2QJE4oERfZJSw1vLUU0+x4IprOey7D/D7huNpmnoiqds2cfXEf+MZZ6C2HA629L6hAy2M3/oA377m6uEJXEQkgm9fezXJWx+M6p5FbTkTaOGUU05h1xFncm/racz8RS3HffFnrK54kNbW1uEJWmSYKdEXcblph2VgjjyFhGOyOOuss3jgjj9gm97npLc3sM4Hu/76E37zg2/T0NDAxRddSPLGG3r+4DzQQvLDN3DBvFzmzJkzvAciItLB3LlzWTDPF9U96+KLLqSxsZGnn36aJ376Bc4OVGF2vkr9lLl86V/jGLPguk5P2xVxCz0wS8SFdgUC3HTnRv74t1d473O/gdSp8NZz8OY34EAzr954IZmZmZ3WMcawZvUqll65jPV3FtIyez5tM88NdmJr3kXCy48xfusDXDAvlzWrV+nBMyLiqIHesz7xsY/y+Mc+ynvvvcf//eEvrFz/DLQeat/ujlkX8Yu//oNjMz/MsUd9yKnDE4kJJfoiLtHU1MSDDz7IXXfdxf2tp9D24U/CYR54bQu8vAn8T7Yv2zXJD0tKSuKONbdRXV3NL371GzasvYqWpkaSU9OYv2AB3/nFvarJF5G4MZh71tSpU/nZ977ByhUdKi0mTofZn+LBpDQevLGG6XvXknfKNFZc+mmOyDhsmI5KJHaU6IuMQNOnT2fHjh2hP04gcfZ5tHk/gv1/X4OWABzzFieMeZ9zj03ht/f/vF/bNsYwd+5c1t3x59gHLiISYzG9Z+3ZDr/PY/ZnLuOltgy2zziLW95O5bcf+zy5man8q+bfNO56H9qCbfozMjL0HBGJa0r0RUaYt99+mx3NFs79Osw4EyZNp/XQAXi9llPP+ChXXJhLfn4+RxxxBJs2beLuO9Z88KUAPcVWRKSjjIyMzvfIaVO49TtLOf3006m45z5+u+ERnnmzjodfbYGzlsJpnwP/U/B6DTveqMNaq6aMErfUGVckRqy1bN68mfwlX2DCRA8JiYlMmOhh0SWXsWXLFqy1A9pu095mbl73CB/91m855pzFHHnkkWAMnPhpeG8bPHQjlOax7VeX8vTjD3DttddyxBFHtK/f8UmSeoqtiEhnPd0jJ02axJcuX0rtulvY8fablJWVwZvPwLbN4D0DPvs/cFUFyZeX8uVly6ioqKChoWHAcQzVZ4iMbqrRF4mBgwcPsvTKZdy/sYp9s+fTtqQMUtJpbt7F3S89yob5n2fBPB9rVq8iKSmpx+10apJz0mcx3rnYI0+D8WlgToC2Z5gwYQJ79+yAWy+Etg86kB177LFDeowiIqPVlClTWLZsGQUFBfDWM4CBw46DY3LYPzaFP/7zj/zxD3+Ai29mUnISZ2Qk8NS6W9njf7Z9G70184nVZ4hIV6rRFxkka23wBr35JZqXlNKWszg4yk1CIqROpS1nMXuXlHLf5hdZeuWybrUy7zXs4vd/rWLh9X9gx1FnfzDjtIXYjBPg1X8wY+v/49qxT7DphktoaGgINr/pkOSrOY6IyND74F5r4d1XyHjjMWpv/iorV67kk2efTcKOl9ndNo5HDhzHngt/CQVr4bQLAdixYwfb3ny72zYH+xki0hvV6IsM0pYtW1i/sYrmJaWQlBx5oaRkWuZdz11/uJQN99/DDTfcwJ3PBXj2YAYtk44KreeFkz8LNWuD61R8F/bt4d1332XatGmdNqfmNyIiw6+ne29WVhbf+973aGpqCvaNeuRh/vTYs3BMNuxvCi40MQPvTc+TFNjAUWOayPnQBOZnebn2ikUE9rXBl/9fn58h6+8spLq6mrlz5w7REYrbqEZfZJBu+vUttMye3/MNOiwpGbIXsaflIN/4xjfY/EI9LQdbSdz6CMf9Zy1LDz4Ct33xg+X37QHoluSLiEh8Sk1N5fzzz+e2m4vhhY2w4afwYlVwZlsrY5+7n4PNjfiTj2dd03F88W8JBGwKZOdH9RnSMns+N918y9AfiLiGEn1xVNfOR+eee96I6Hw0ffr09qcolq9bS9sJ50W34qzzoK2VK6+8kt9f/lFqv3MWzRt+zivrS7n95hu7NcFRkxwRkZGp2/18QiItlb/hxRs+y21zGlhy8DFmvHoPNO6E2blRbbNt5rmsu+uuuH+K70j9bHcj1zTdMcZ4gBVAfWhSprW2aLjWl/4brs5HnTq4Erz5btu2jebmZvbu3Utzc3P7/xcuXMiuXbval01JSeGiiy5i53vv8U6gmXf3JRB461Vadu6AKccGh7is+G7wSYzRSEmH1oP88Y9/jDhbTXJERNyhp/v5CSecwAknnMDlS4N/G3Nr/z5D2lphXhHsfocdgXc48uijmTZlCtOmTWPq1KlMmzaN1atX09TU1L5aeno6jzzyCCkpKaSkpDBhwgRSUlJITk7miCOO6PYZOZjPInUsji+uSfSBcqDQWusHMMZ4jTGV1troviYPfn3ph66djzr9ZBnufHTqBdy38QaWXrmMO9bc1mmc4oyMDN59990PVklN5Vvf+hZ79uyhsbGRxsbG9v93vIFBsENUSkpK54BMAqR4YMJUmDgu+NCUlHSafd/kz0yHzMNhbGidx34DO9+E1gOQmARJ46F5V7DzVF+ad8EY3dhERCRkTFL/PkPGJsORp8IsHxxo5u1bF/L2m29C7rdg3Mnw5jsw9wrYvR0aXodtW9i1axdz5swF+q5J37FjB5/+9KdJS0tj4sSJpKWltf//Jz/5CY2Nje3Ldv1SMNjPdok9VyT6xpg8wB9O0gGstf7QT1t51tqKoVxf+q8/HVjvu2MZN954I/v37+e5557j+eef75TkAzQ1NXHDDTd038a4VPAcAeMnBmtCJkyBPdsZv+N5klMn0nT+z2hLSad13MRgsg9QvRb+XgYHWsDzIU6dkcHRk8Zw3NSxzDx8Ep/85ipmHXUvBP4Ld10DY8bB1kqYu6TP4054+THyFi/uZ2mJiIhb5ecv4u6XHg2OttOXrZXQehD+cEmwoiklnddee4333nuP2+reZ8v2g/x371Tebjox+Pn33jbYtgWA1CvLaEmZRkLzLtjbQFvjTlrffhGeuS+47ane4Lb37eGRyiqwbX2Gs2PHDubPn89JJ53ESSedhLWW+zdW0rykTB2L44QrEn1gMVAZYXolUAj0lagPdn3pp/50YG2Z9Vn+58clMCkDxk6AQx1qxE9bCJOCifyHT8nmQMI4MlMPsey4A6SlpXHJ35PZc7DLNl96lJbXqgG46E/P4UkewxETx3H4xLF87YsXw/vbgssd2gdrvsTTEdoSdnqS4qH9UFcBp1/Y+/EcaGH81gf49s/v6f2YRURk1Pj2tVez4fyL2HvqBX1+hvDvu+HQgeDfrQfJSDEcc8wxHHPMMWRnf7CoMSb4/JXxae3TVl4xj+e3N/Hf3Qd4p3E/7+w5wGlHXMyDy0KJ/sIbYNLh7cunJLaRk9bE59JeY8+ePdyz+yie/ncttOyB/XvhQDO8/xobNmxgw4YNMPmY4C8OOYv61bF47V/W9LfIpB/ckuj7gNII0/1AzjCs7xrWWlrbLIkJwY4+zQdaadrfysE2y6HWttC/lhMyJgBQ/14z/91zgP2H2tpfAJ8/5TAAHtz6Hlt37OVAq2XfwTaaD7aSkpTIgw+sD7bbi8asXKhZB5f8FoBTpo3h2RWfDM475XyYeDjs20Pq9GOYnJLEJ7weLv7UDAB+Pe0drr326+zZ8SY074a973NY6ljgpwDcffnJnXZ1Q9Or7Ah0bqsYSdefKi9ZegX3bbyBlnnXR77BHWgh+eEbuGBeLnPmzInuuEVExPXmzp3Lgnm+qD5DFn5+YVTNXdoro/Y1tv999ceO7LactZbDvx9a9pFfwITJpE37EN/+3x/xfvNBTpw+gcKzLgTgrpKnYNYUGJ/6wQaeWc/dt/6UZ559jhsaPx58kGM/Ohbfd2chAI37DvGrv7/JuMQExo0JvwxnHjOJWRkTaNp/iH+9trt9XlKCISkxgaM840hPSWL/oTbe23uQMQmGpERDUoJhTKJhXGICCQmju2nQiE/0Q51oPUCk504HQvOGbH2A/+6Fz90WfPpduPL322cfzce9Hp79bxPff8iPxWLtB63jfvCpGcw9eiJPvb6bH2zc1m3+/13wYU45IpXKlxv4cdU2rIU2G3xTtlm4fcksZh42gbuffZefPfo6rW2WttC8Nmt56MuncszkZMqefJuVj71Om7W0tkFrm6XVWp7/zhkcljaWH1du42dVr9MaSvDbQgHsXflJUsYm8j8b6vn139/qUmbQ9otzAfjZo6+zess7neanjUtsT/TX1Gxn3TPBZjYJBlLGJjJj8nha9jb2q/OROdDMA18+hbRxYzgsNYlP/ip0Y/pzIdg2MjIyqFvVvfPQ5XMP5/Ino/9BZiAdkIwxrFm9iqVXLmP9nYW0zJ5P28xzg8fXvIuElx9j/NYHuGBeLmtWr1J7RBERaTcUnyHRfpYZY6JeduvyM4ODW+x8L9hnLSmZqZM9fP7zv+SChRdy4nM7WfzLvf36bD/QHOwwvHPvQa7fuK3bIjdf+GFmZUzA//4+Pl32TLf5qxefwBVzj6D2zUY+ekttt/kVXzyJi045jMqXG5j/x2dINIbEBENiAoxJMFR88WTOOS6dh158n6vufokEY0gwkGAMicZw5xdO5PQj07jv+Z384OFt7fMSTLDs7rh0NplTU6h45l1+/fc3MQYM4flwx6UnMn3iOO6o286fqt/BYILLhJa767ITmTh+DLdXv0PFs8FcyWBC5wbu/uJJjElMYNVTb/PQS5HS1L6ZkT7EkTHGS3CknGxrbV2XeT6CzW/SrbWBWK5vjCkACgASD8vMPuYrZYTODQYomGWYe5jhpYDl/56x7TNMaP5XTzScMsXw7PuW32+1oZP+wfrXnGz48CRD7U7LX175YH5CaIFrTzIcmWrYvMNy3+sWAyS2Xzxw9UmGqeMNT+6wPP5fG5xHMNlONLBsliE1ybDlXUvde8H54XkJxnBxJoxNDMZXvwfGJATnjTHB/59zRPAi37bH0rAfkhI+eI1NhKNTg0fTcih47EkJ4fiC0+fNX8D+L/wxus5HjTsZf0cBDz1wX9/LOshay0svvcTav97H5iefZH9LE+OSUznrI2ex6PMXcsIJJwx7TE1NTaSmpva9oERF5RlbKs/YUVnGlhPlGY+fIf017/wL2H/pHwb02d7aZjnQBgc7vFKTIDXJ0HLI8uqe0PRWOGThUBvM9MD0FEPDPsu/dnwwvdVCaxt8bDocnWZ4s8my8U1La6jSNPxaeKzhmLRgrnbfa/aDeQQrbq+cGcy1anda7tlmsQS3DcHlvnWKYXqK4e/vWO557YMK2/C/P8w2TB5v2PimZf3rnSt0AX5+ZjAX++u2YHzhmeFlfvdxw5gEw19esTz+dnD/2350Xq21NurWJkr0B7k+wMyZM+3LL788mMMYNd577z1++ctfUvJ/v+ZQzpLoOrDWrCXvqBa14xuATZs2cfbZZzsdhmuoPGNL5Rk7KsvYUnkOzKJLLuPut1Oi61i8+Q5mvPtPKh96gMzMzKEPziWMMf1K9Ed80x2Jf13HsQ9LePoe2tSBVURExBX61bG4roJtLbs57rjjOs0a7Dj+0pkbnozbZ6Ol3mrjY7C+9CFSkv/Pf/6TRReeT/LGG+BgS+QV1YFVRERkxAh3LI7ms/38T53H5Zdf3m12pJxBBm7EJ/qhJDwAeCPM9obmDdn60rt9+/ZFnP6Rj3yENatXsfCMWUy4s5CEmrXBx4C3HoLGnSTUrCXlrgIWnjFLHVhFRERGgHDH4mg+2/9afhe33Xab0yG7nlua7tQAkyNMzwSqhmF9iaCtrY0rr7yyx/lJSUncseY2qqur+cWvfsOGtVfR3LSHlNSJzF+wgO/84l7V5IuIiIwgsfhsV+fy2HFLol8O5AJdB2b3ASuHYX2J4Prrr+fOO+/sNr3j2PTGGObOncu6O/4MqAOUiIjISNefz/ZOD6AMueSSS7jnnntITEwc6lBdb8Q33QGw1pYB3tAIOgAYY7KABmttp0HUjTH1xpjaga4v0Vm9ejU//elPSUxM5KGHHsJa2/5SJxsRERGB4Jj/4fzgpZdeIj09nfXr1/PNb37T6dBcwS01+gDnASuMMfWhvzOttZEezxYgcgfcaNeXPlRVVVFYGHza3a233sq8efMcjkhERETi3cyZM7n33nvJzc3lN7/5DZmZmVx77bVOhzWiuSbRD3WqLYpiuezBrC+9e+GFF7jooos4dOgQ3/3ud9sTfhEREZG+fOITn2D16tV84Qtf4Jvf/CbHHnssCxcudDqsEcsVTXckPmzfvp3Pfvaz7Nmzh7y8PG688UanQxIREZER5tJLL+WGG27AWssll1xCTU2N0yGNWEr0JSb27t3LggULeOONNzjzzDNZs2YNCQm6vERERKT/vv/973P55ZfT3NzM+eefz+uvv+50SCOSMjEZtNbWVi699FJqamqYMWMG9913H8nJvTwRT0RERKQXxhhKS0s599xz2bFjB/Pnz2f37t1OhzXiKNGXQfvOd77DfffdR3p6Ohs2bOCwww5zOiQREREZ4caOHcvdd9/NrFmzeOGFF8jLy+PgwYNOhzWiuKYzrgy/6dOndxv79oQTTnAoGhEREXEbj8fDhg0bOOOMM6iqqmLs2LHt8zIyMjRkdx9Uoy8D1jXJ37Vrl0ORiIiIiFsde+yxPPDAA92md81DpDsl+iIiIiIS1+bMmeN0CCOSEn0ZkFdffdXpEERERESkF0r0ZUB+97vfdZuWkZHhQCQiIiIyGnTNM5R39E2JvvRbc3Mzq1evBqC6uhprLdZadYgRERGRIbN9+3YCgQATJkwA4LHHHnM4ovinRF/67Y477iAQCHDGGWeQk5PjdDgiIiIySkyaNInLLrsMgN/+9rcORxP/lOhLv1hrueWWWwD42te+5nA0IiIiMtqE84/bb7+dPXv2OBxNfFOiL/3yr3/9i2eeeYZp06aRn5/vdDgiIiIyypx00kl88pOfpKmpiT//+c9OhxPXlOhLv4Rr87/85S8zfvx4h6MRERGR0Shcq3/rrbdirXU4mvilRF+i9s4771BRUUFCQgJXXXWV0+GIiIjIKHXhhRdyxBFH8OKLL/L44487HU7cUqIvUVu1ahWHDh3iggsu4Oijj3Y6HBERERmlkpKSKCwsBD5obSDdKdGXqBw8eJDS0lIArr76aoejERERkdFu2bJljBkzhvvuu4833njD6XDikhJ9icq9997Lf//7X0444QTOPfdcp8MRERGRUe7www8nLy+Ptra29spI6UyJvkTl1ltvBYKdX4wxDkcjIiIi8kGn3FWrVrF//36Ho4k/SvSlT88//zxPPPEEqampLF261OlwRERERAD46Ec/yqmnnsrOnTupqKhwOpy4o0Rf+hSuzV+6dCkTJ050OBoRERGRIGNMe62+OuV2p0RferV79+72h1F89atfdTgaERERkc4uueQSPB4PTz31FLW1tU6HE1eU6Euvbr/9dvbu3cs555zDiSee6HQ4IiIiIp1MmDCBK664AvigFYIEKdGXHrW1tXXqhCsiIiISj8KtDu68807ef/99h6OJH0r0pUePPvoo//nPfzjyyCNZuHCh0+GIiIiIRHTccccxb9489u3bx+rVq50OJ24o0ZcehWvzCwsLGTNmjMPRiIiIiPQs3Prgd7/7Ha2trQ5HEx+U6EtEr7/+OuvXrycpKYlly5Y5HY6IiIhIrz7zmc8wY8YMtm3bxkMPPeR0OHFBib5E9Pvf/562tjby8/PJyMhwOhwRERGRXiUmJvKVr3wFUKfcMCX60s2+fftYtWoVAFdffbXD0YiIiIhE58orr2T8+PFs3LiRV155xelwHKdEXzqZPn06ycnJ7T3WP/e5zzkckYiIiEh0pkyZ0v7/448/HmMM06dPdzAiZynRl0527NjR698iIiIi8Wzfvn2d/h7NuYwSfRERERERF1KiLyIiIiLiQkr0pZOUlJROf2vEHRERERlJuuYuozmXUaIvnZx44okAVFVVYa1l+/btDkckIiIiEr3t27ezc+dOAJKTk3nzzTcdjsg5SvSlXUtLC//+979JSEhg7ty5TocjIiIiMiBTp07lwx/+MC0tLTzzzDNOh+MYJfrSrra2lkOHDnHyySeTlpbmdDgiIiIiA/aRj3wEgCeffNLhSJyjRF/a/etf/wLgrLPOcjgSERERkcEJ5zPh/GY0GuN0ALFgjPEAK4D60KRMa21RP9avB4qBdaFJPiAXKLLWBmIXaXwLf+NVoi8iIiIjXTifGc01+q5I9IFyoNBa6wcwxniNMZXW2two1/cCpaEXgB/IHU1JvrVWib6IiIi4xoknnkhaWhqvv/4677zzDocffrjTIQ27Ed90xxiTB/jDST5Ah4Q/L8rNFAHZBGvxs621mR23Nxq89tpr7Nixg6lTp3Lcccc5HY6IiIjIoCQmJnLGGWcAo7dWf8Qn+sBioDbC9EqgMNqNWGvrrLVV1tq6mEU2gnRsn2+McTgaERERkcEb7e303ZDo+wg2tenKD+QMcywjlprtiIiIiNuM9nb6IzrRD3XC9QANEWYHQvP6s70sY4xvsHGNREr0RURExG3OPPNMIDiE+IEDBxyOZvgZa63TMQyYMcZLcKSd7K5NbkIJeyWQ3lenWmNMObAWqLPW+o0xxUC9tbasl3UKgAKAadOmZa9bt66nReNeS0sL559/PgAPPPAAycnJjsXS1NREamqqY/t3G5VnbKk8Y0vlGTsqy9hSecaW0+X5xS9+kTfeeINbb72V2bNnOxZHLJxzzjm11tqoW6y4ZdSdwSq11laF/7DWFhlj6o0x/o7TOwp9CSgDmDlzpj377LOHJ9IhsGnTJtra2sjKyuIzn/mM47GM5LKMNyrP2FJ5xpbKM3ZUlrGl8owtp8vT5/OxevVq9u/fP+rO64huuhMrPSTzFQTH1nc9NdsRERERtxrN7fQdr9EPPayqPwLW2uzQ/yO1ze9kEGPh1wPLB7juiKJEX0RERNxKib6DrLWZg1g3YIwJEHzgVddhMb0EO+T2KvRFo8haWzHQOEYyPShLRERE3GzWrFlMmjSJt956i7feeosjjzzS6ZCGjRua7tQAkyNMzwQitq/vIkD3Lwnh9V0/pn59fT3vvfceGRkZzJgxw+lwRERERGIqISGhffSd0Var74ZEv5zgE2278hEcSacva3t4Cm4eUDqYwEYCPShLRERE3G60PjhrxCf6odFvvKGhNoHgePhAQ9fmOKGRdLo+RbcqNFRmx+WWA/7ehtd0CzXbEREREbcbre30HW+jHyPnASs6dOzNtNZGquUP0KUDr7W2zhhDaOx8CD5kq76H9V1Hib6IiIi43RlnnIExhrq6Ovbt28f48eOdDmlYuCLRD42sUxTFctk9TK9jFLTH76qxsZHnnnuOMWPGkJMT9bMXREREREaUSZMmceKJJ/L8889TV1fHRz7yEadDGhYjvumODNyWLVtoa2vj9NNPd/RpuCIiIiJDbTS201eiP4qp2Y6IiIiMFqOxnb4S/VFMib6IiIiMFh0TfWutw9EMDyX6o1RbWxtPPfUUoERfRERE3O/4448nPT2dd955hzfeeMPpcIaFEv1R6j//+Q8NDQ0cccQRHH300U6HIyIiIjKkEhISRl07fSX6o1THZjt6UJaIiIiMBqOtnb4S/VFK7fNFRERktFGiL6OCEn0REREZbebOnUtCQgJPP/00zc3NTocz5JToj0K7d+/mhRdeYOzYsWRlZTkdjoiIiMiwSEtL4+STT+bQoUPU1NQ4Hc6QU6I/Cm3evBlrLVlZWaPmEdAiIiIiMLqa7yjRH4XUbEdERERGKyX64mpK9EVERGS0Gk0PzlKiP8roQVkiIiIymh133HFMnTqVd999F7/f73Q4Q0qJ/ijz4osvsnv3bo466iiOPPJIp8MRERERGVbGmFHTfEeJ/iijZjsiIiIy2inRF1dSoi8iIiKjnRJ9cSUl+iIiIjLazZkzh8TERJ599lmampqcDmfIKNEfRRoaGnjxxRcZN24cp59+utPhiIiIiDhiwoQJnHrqqbS2tlJdXe10OENGif4osnnzZgBycnIYO3asw9GIiIiIOGc0NN9Roj+KqNmOiIiISJASfXEVJfoiIiIiQaPhwVlK9EeJ1tZWPShLREREJGTGjBlkZGTw/vvv88orrzgdzpBQoj9KvPDCCzQ1NXHsscdy+OGHOx2OiIiIiKNGw4OzlOiPAtOnT+fUU08F4LXXXmP69OkORyQiIiLivKqqKgAuv/xyjDGuy5GU6I8CO3bs6PVvERERkdGo6xj6bsuRlOiLiIiIiLiQEn0RERERERdSoj8KHHbYYZ3+zsjIcCgSERERkfjRNSdyW46kRH8UeOKJJwDIzMzEWsv27dsdjkhERETEedu3b+eiiy4C4I477nBdjqREfxTYtm0bAMcee6yzgYiIiIjEmXB+9Nprrzkax1BQoj8KhC/cGTNmOBuIiIiISJwJ50fhilE3UaI/CoQvXCX6IiIiIp0p0ZcRLVyjr6Y7IiIiIp2p6Y6MaKrRFxEREYksnOi//vrrtLa2OhtMjCnRHwXUGVdEREQkspSUFA477DAOHjzIO++843Q4MaVE3+UaGxt5//33GT9+PNOnT3c6HBEREZG449Z2+kr0XS7c3uyYY47BGONsMCIiIiJxyK3t9JXou5yG1hQRERHpnVtr9Mc4HUAsGWPygMXW2vx+rucBVgD1oUmZ1tqiGIfnCHXEFREREeldOE9yW42+KxJ9Y0xp6L9eYPIANlEOFFpr/aHteY0xldba3FjF6BR1xBURERHpXThPcluNviua7lhrC621hQQT9n4J/QrgDyf5oe35O8wb0dR0R0RERKR3bq3Rd0WiP0iLgdoI0yuBwmGOJeZUoy8iIiLSu6OPPhpjDG+++SaHDh1yOpyYUaIPPsAfYbofyBnmWGJONfoiIiIivRs3bhxHHHEEra2tvPXWW06HEzOjOtEPdcL1AA0RZgdC80asXbt2sXv3biZMmMCUKVOcDkdEREQkbrmxnf6oTvSJouNu6MvAiNSxNl9j6IuIiIj0zI1DbBprrdMxxIwxpoDg6DnZUS7vJTikZra1tq7LPB/Bdvrp1tpAD/sqCP15EvD8IEKXD0wF3nM6CBdRecaWyjO2VJ6xo7KMLZVnbKk8Y2emtTYt2oVdMbymE6y1ZUAZgDGmxlo74tvzxwOVZWypPGNL5RlbKs/YUVnGlsoztlSesWOMqenP8o4n+saY+r6X6iQQbY19FCK1ze8kUm2+iIiIiEi8czzRt9ZmOrjvgDEmQPBBW3VdZnsJdsgVERERERlxRntnXIAaInfKzQSqotxGWezCGfVUlrGl8owtlWdsqTxjR2UZWyrP2FJ5xk6/ynJUd8btsE6utTa/y/RaYKW1tiLGYYqIiIiIDDm31eh76GXse2NMfSiBbxfqVOsNjcATXi4LaFCSLyIiIiIjleNt9GPBGFNMMMFfBHiMMeUEO9qWdhk2M0DkDrjnASs6dAzOtNbmDl3EIiIiIiJDy1VNd+JVqHmQ31obbZt/6SD0a0thh0leoMha63copBEt9IuVD5gCZAF+guUZcDKukcwYkwcs7toEULoLPYRwBcFnmECwYqXIuYhGNl17saN7Y2zps3vo9CevdEWNfjwLfaiVAvqFYABCN4q8jolA6IOt1hiTrRtG/4Q+yHKstSUdpi0HdhljMlWe/WOMKQ3910sUT9oWAMoJ9qXyQ/A9boyp1K+o/aNrL7Z0b4wtfXYPnf7mlW5rox+PFqFhOgejsOuEUN8JT6R50idfqF9Ku9AHWx3BG4f0g7W20FpbSDB5lT6EPuj9HT/kOyT8eY4FNgLp2os53RtjS5/dQ6dfeaUS/SFkjPER/RCd0rPFEaYF6KXjtfSoMHRddlVF8CdrkaG0GKiNML0SffiLs3RvjD19dsfYQPJKJfpDy6ufpwbHWlvUdbjU0M9WHlSTNVDevhcRGRI+gu2eu/IDOcMci0hXujfGiD67h0y/80q10R8ixpiCrj8DSsysAsrUubn/enkSdaSnQ4vETIcP+UgjnwVQLZ84SPfGYaHP7kEYaF6pRH8IhDqhqCY/hjqMhjAHWKtnHMROKAHLAzRqhwylPjuMGmM8GuFE4oXujYOnz+7YGExeqUR/aHTr1CODE3oeQl3oYi8yxkxWGcdMMVChG7CISCe6Nw6SPrtjZsB5pRL9GAuNHLHO6TjcKtQ2rdAYsytU+1fS50rSo1DHnpyubSlFREYz3RtjS5/dAzfYvFKJfhcdno4brUD4RhD6mQ/99PyBwZRnH8oI1raMqptFLMszdL0WEXwy9KgzhNemRBapbX4nundKPBjt98YhNio/uwcqFnmlEv0ueumQE40VgMcYM6fLdA/Bn6xygcrR1BFlMOUZusDLCT5Jr2tnqPdDy4yqkY0GeX12tQrIH63JVYzLUvpgrQ0YYwJE7tzoRc8bkfgxqu+Ng6XP7pgadF6pRD+GenqMe+jpesWjKcGPES/BTjyREoMpoX/7rCWU7owxxXR5tLsxJivCTVkklmqI3Ck3Ez1zROKA7o0xoc/uGIlFXqlx9CVuhW6sJT10hPIBVapx6T9jTAHB0Q+61qZoHHMZauVEfmy7D1g7zLGIdKJ7Y2zoszu+KNEfPh6nAxihqkM333ahTlJe9CTNfguVXWbo/1mhly80Xe3PB86D3uN9Co0a4Q2NwAG0D7/XoJFNBsyDrr1B070x5vTZPfQ80SxkrLVDHMfoFfoJ0EtwHF4/wZ+mi9UurX9CicBigm37phAs02WqEeg/Y0xvb/gya61uwP0Qeo97gEWhfysI/iRdqp/6Iwu1310BhDtDZ/b087T0TNdebOneGHv67I69geSVSvRFRERERFxITXdERERERFxIib6IiIiIiAsp0RcRERERcSEl+iIiIiIiLqREX0RERETEhZToi4iIiIi4kBJ9EREREREXUqIvIiIiIuJCSvRFRERERFxIib6IiIiIiAsp0RcRERERcSEl+iIiIiIiLjTG6QBERGR0MMYsBwJALlAENAAFodlTrLVFDoUmIuJKSvRFRGTIGWOKgZXW2oAxxg+UA1XW2iJjjA+oNMastdbWORupiIh7qOmOiIgMqXAib60NhCb5gSygPvS3ByhTki8iElvGWut0DCIi4mLGmDxrbUWHvwuAUmutcTAsERHXU6IvIiLDyhhTDnittdlOxyIi4mZquiMiIsPNB6x1OggREbdToi8iIsPGGOMl2Ca/qsv0LEcCEhFxMSX6IiIypIwxpaEEH6AQoGPH21Cb/YADoYmIuJoSfRERGTLGmDyCY+V7QpPeD033dvjXY631OxKgiIiLqTOuiIgMGWOMBygGagkm9CWh4TbzQ9Ow1pY5F6GIiHsp0RcRERERcSE13RERERERcSEl+iIiIiIiLqREX0RERETEhZToi4iIiIi4kBJ9EREREREXUqIvIiIiIuJCSvRFRERERFxIib6IiIiIiAsp0RcRERERcSEl+iIiIiIiLvT/AfrAPVpM9R/MAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Set the target interval and number of (equidistant)\n", "# points we wish to sample\n", "xStart = -3.\n", "xEnd = 3.\n", "nPoints = 11\n", "\n", "# Create the list of points based on the specifications \n", "# above (uncomment example function 1/2)\n", "x = np.linspace(xStart, xEnd, num=nPoints)\n", "# Example function 1\n", "y = np.exp(-pow(x,2.)/2.)\n", "# Example function 2\n", "# y = 1/(1+pow(x,2.))\n", "\n", "# Generate a finer dataset in order to assist plotting\n", "x_fine = np.linspace(xStart-1, xEnd+1, num=64)\n", "# Example function 1\n", "y_fine = np.exp(-pow(x_fine,2.)/2.)\n", "# Example function 2\n", "# y_fine = 1/(1+pow(x_fine,2.))\n", "\n", "# Method 1: Home-made Lagrange interpolant\n", "def home_lagrange(xs,ys):\n", " n = len(xs)\n", " p = np.poly1d(0.0)\n", " for j in range(n):\n", " current = np.poly1d(ys[j])\n", " for k in range(n):\n", " if k != j:\n", " current = current*np.poly1d([1,-xs[k]])/(xs[j]-xs[k])\n", " p = p + current\n", " return p\n", "\n", "home_p = home_lagrange(x,y)\n", "\n", "# Method 2: In-built functionality\n", "inbuilt_p = lagrange(x, y)\n", "\n", "\n", "# Plot results\n", "plt.rcParams.update({'font.size': 22})\n", "plt.figure(figsize=(12, 8))\n", "\n", "plt.plot(x_fine, home_p(x_fine), \"k-\", linewidth=2.0, label=r\"Home-made Lagrange interpolant\")\n", "plt.plot(x_fine, inbuilt_p(x_fine), \"ks\", markersize = 4, label=r\"In-built Lagrange interpolant\")\n", "plt.plot(x_fine, y_fine, \"--\", color = (0, 0.447, 0.741, 1), label=r\"Target function\")\n", "plt.plot(x, y, \"ko\", markersize = 12, markerfacecolor=(0, 0.447, 0.741, 1), label=r\"Initial data\")\n", "\n", "plt.legend(loc=\"upper left\")\n", "\n", "plt.xlabel(r\"$x$\")\n", "plt.ylabel(r\"$y$\")\n", "plt.xlim([xStart-1, xEnd+1])\n", "plt.ylim([-1, 3])\n", "plt.grid()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Exercise 1**: modify and extend the above code such that you reproduce the results from Figure 3.4 in the notes.\n", "\n", "**Exercise 2**: practice finding the coefficients for simple functions. There are a number of good [resources](https://pythonnumericalmethods.berkeley.edu/notebooks/chapter17.04-Lagrange-Polynomial-Interpolation.html) out there." ] } ], "metadata": { "celltoolbar": "Slideshow", "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.8.5" } }, "nbformat": 4, "nbformat_minor": 2 }